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A Random Search Algorithm for Laboratory Computers 


One of the major accomplishments during this reporting 
period was the completion of the year's second doctoral thesis. 
Dr. Arye R. Ephrath completed his degree requirements in June 
of 1975 and his thesis j abstracted here, has been made into a 
NASA Contractor's Report #GR-137759, 


PILOT PERFORMANCE IN ZERO-VISIBILITY PRECISION 

APPROACH 

Arye R. Ephrath 

PH.D. Thesis, Department of Aeronautics and Astro- 
nautics, Massachusetts Institute of Technology, June 

1975 

ABSTRACT 


To gain a better understanding of the pilot's short 
term decisions regarding performance assessment and 
failure monitoring, the performance of fifteen air 
line pilots who flew simulated zero -visibility landing 
approaches is reported. The approaches were flotra 
with different degrees of automation an!| at different 
workload levels, as induced by simulated wind distur- 
bances and measured by a sensitive, non-loading 
subsidiary task. 

The results indicate that the pilot's mode of parti- 
cipation in the control task has a strong effect on 
his workload, the induced workload being lowest when 
the pilot acts as a monitor during a coupled approach 
and highest when the pilot is an active element in 
the control loop. In addition, a very marked increase 
in workload at altitudes below appro aximately 500 feet 
A6L is documented at all participation modes; this 
increase is inversely related to distance-to-go - 

The effects of workload and participation mode on 
failure detection are separated. The participation 
mode is shown to have a dominant effect on failure 


detection performance, with a failure in a 
monitored (coupled) axis being detected faster 
than a comparable failure in a manually con- 
trolled axis. 

Touchdown performance is also documented and the 
findings of previous investigators are supported, 
namely that the conventional instrument panel and 
its associated displays are quire inadequate for 
zero-visibility operations in the final phases of 
the landing approach. 

A paper summarizing some of the work presented in this 
thesis was delivered by Dr. Ephrath at the Eleventh Annual 
Conference on Manual Control. The paper, entitled "Detection 
of System Failures in Multi-Axes Tasks”, is contained in 
Appendix A. 

A research note in preparation by Dr. Curry and Dr. Ephrath 
is contained in Appendix B. Its title is "Changes in Pilot Work- 
load During an Instrument Landing" 

A paper on modelling pilot performance in failure detection 
tasks was presented at the Eleventh Annual Conference on Manual 
Control by Drsv Curry and Gai. The paper has since been revised 
and submitted for publication. The abstract of the paper which 
was presented is given below and the text of the paper is given 
in Appendix G. 


FAILURE DETECTION BY PILOTS DURING AUTOmTIC LANDING: 
MODELS AND EXPERIMENTS 

Eli G, Gai and R,E. Curry 

ABSTRACT 


A model of the pilot as a monitor of instrument 
failures during automatic landing is proposed. 

The failure detection model of the pilot consists 
of two stages: a linear estimator (Kalman Filter) 
and a decision mechanism which is based on sequen- 
tial analysis. The filter equations are derived 
from a simplified version of the linearized dynamics 
of the aircraft and the control loop. The percep- 
tual observation noise is modelled to include the 
effects of the partition of attention among the 
several instruments. The final result is a simple 
model consisting o£ a high pass filter to produce 
the observation residuals j and a decision function 
which is a pure integration of the residuals minus 
a bias terra. 

The dynamics of a Boeing 707 were used to simulate 
the fully coupled final approach in a fixed base 
simulatoi which also included failures in the air- 
speed, glideslope, and localizer indicators. Sub- 
jects monitored the approaches and detected the 
failures; their performance was compared with the 
predictions of the model with good agreement between 
the experimental data and the model. 


During this reporting period, the data collected during 
the principal investigator’s year ac Ames Research Center was 
analyzed. A discussion of the experiments performed, the data 
analysis, and the results of these experiments are contained 
in a paper presented at the Eleventh Annual Conference on Manual 
Control. The paper is abstracted here and presented in full in 
Appendix D. 


EXPERIMENTS IN PILOT DECISION-lfeKING DTJRING SIMDLATED 
LOW VISIBILITY APPROACHES 


R.E, Curry, J.K. Lauber, & G.E, Billings 
ABSTRACT 


Despite a vast accumulation of operational experience 
with the conduct of low visibility instrument approaches, 
little is understood about the decision making behavior 
of pilots who fly these approaches. Likewise, there 
is little information regarding the man, system and 
task related factors which influence this decision- 
making behavior. Such information is essential for 
the rational design of new systems, or for the redesign 
of existing systems in order to correct known deficiencies 

A major problem which has inhibited the study of pilot 
decision-making behavior in the laboratory has been the 
unavailability of tasks which incorporate the essential 
cognitive features of the real task, and which include 
those motivational or stress-inducing features known to 
influence decision-making performance. This paper 
describes a task which was designed to simulate the 
cognitive features of low visibility instrument approaches 
and to produce controlled amounts of subjective stress 
in pilots serving as subjects in experiments using the 
task. 

Pilot behavior during low visibility instrument approaches 
can be analyzed into at least two major categories: one 
is the continuous closed-loop manual tracking behavior 
necessary to control the aircraft, and the other is the 
cognitive, decision-making behavior required to make the 
decision to continue the approach and landing or to 
execute a missed approach. It is the second category of 
behavior with which we are concerned here. 

The difficulty of a decision making task is, in part, 
determined by the uncertainty of the data used to make 
the decision. For example, the decision to "go-around" 
is a relatively easy one if, at the missed approach 
point, there is nothing to be seen outside the aircraft, 
or if the approach lights and runway have been clearly 
visible for the last two miles of the approach. It is 
when the approach lights or runway are barely visible, 
and then only intermittently, that the decision becomes 
more difficult. 


A second class of variables which, are known to influence 
the outcome of decision-making tasks is best illustrated 
by the various kinds of psychological stressors acting 
upon the pilot. Of particular interest here are the 
pressures perceived by the pilot to complete the approach, 
to make an on-time arrival, to save fuel, and even to 
save "face". 

We have assumed that it is necessary to use a simulation 
task which incorporates both kinds of variables, infor- 
mational and psychological, to successfully study pilot 
decision-making behavior in the laboratory. The task 
discussed below T-?as designed to meet those requirements. 
This paper describes our preliminary experiments in the 
measurement of decisions and the inducement of stress 
in simulated low visibility approaches. 


PUBLICATIONS 


Two major publications resulting from this research 
appeared during the reporting period. Both appeared in 
Behavior Research Methods and Instrumentation and they are 
abstracted below. Pull, copies of the papers appear iu the 
Appendices E and E. A more detailed description of MDNOML 
containing program listings and subroutines is available 
from the principal Investigator. 


A rOLTraOMIAL MAXIMUM LIKELIHOOD PROGRAM (MDHOML) 

R.E. Curry 

Behavior Research Methods and Instrumentation, 
Volume 7(3), pages 305-307, 1975. 

ABSTRACT 


In our research on modelling sensory and decision 
phenomena, we were soon confronted with the task 
of evaluating both old and new models using both 
old and new data. Rather than design an ad-hoc 
estimation program for each new model, as is. typ- 
ically done, we developed an "executive" program 
that provides a general method for estimating 
parameters and simultaneously provides flexibility 
for accommodating new models with a minimum amount 
of programming. Our experience with canned computer 
programs has been equivocal, so we decided to 
provide only the general framework and let the 
user accomplish the objective of estimating the 
parameters for his particular model by writing a 
new subroutine within the constraints of the 
executive program. 

The most common class of distributions for which 
parameters must be extracted are multinomial 
distributions resulting from a stimulus-response 
classification, e.g, binary responses (YES-HO or 
two alternative forced choice methods) , the method 
of successive categories (rating scales), or 
transition probabilities in a Markov chain- Although 
a number of methods exist for estimating such 
parameters, we have chosen the maximum likelihood 
method and have implemented the scoring of Rao to 
adjust the parameters from one iteration to the next. 
We have chosen the maximum likelihood (ML) method 
because (1) it is a member of the class of consistent 
asymptotically normal estimators (CAN) under suitable 
conditions; (2) it will easily handle situations in 


7 



which all the respotties fall into one category; 

(3) there are many situations in which the ma^cimum 
likelihood estimator can be shown to yield unique 
estimates for parameters (other estimation techniques 
may have this property as well) ; and (4) it is the 
only one exhibiting first order efficiency. 

Other authors have successfully used the maximum 
likelihood method with Rao's scoring method in 
signal detection type models. This program differs 
primarily in its ability to estimate parameters 
for a wide variety of models. 


A RANDOM SEARCH ALGORITHM FOR LABORATORY COMPUTERS 
Renwick E. Curry 

Behavior Research Methods and Instrumentation, 
Volume 7(4), pages 369-376, 1975. 

ABSTRACT 


The small laboratory computer is ideal for experimental 
control and data acquisition. Postexperimental data 
processing is many times performed on large computers 
because of the availability of sophisticated programs, 
but costs and data compatability are negative factors. 
Parameter optimization, which subsumes curve fitting, 
model fitting, parameter estimation, least squares, etc., 
can be accomplished on tha small computer and offers 
ease of progamm’ng, data compatability, and low cost 
as attractive features. A previously proposed random 
search algorithm ("random creep") was found to be very 
slow in convergence. We present a new method (the 
"random leap" algorithm) which starts in a global search 
mode and automatically adjusts step size to speed con- 
vergence. A FORTRAN executive program for the random 
leap algorithm is presented which calls a user 
supplied function subroutine. An example of a 
function subroutine is give which calculates 
maximum likelihood estimates of receiver oper- 
ating characteristic parameters from binary 
response data. Other applications in parameter 
estimation, generalized least squares, and 
matrix inversion are discussed. 




DETECTION OF SYSTEM FAILURES IN MULtl^AXES TASKS 


by Arye R* Ephratb 
Hau~Vehicla Laboratory 
Massachusetts Institute of Technology 


SUMtlARY 


The investigation has examined the effects of the pilot’s participation 
mode in the control task on his workload level and failure-detection perfor- 
mance during a low-visibility landing approach. We found that the participa- 
tion mode had a strong effect on the pilot’s workload, the induced workload 
being lowest when the pilot acted as a monitoring element during a coupled 
approach and highest when the pilot was an active element in the control loop. 

The effects of workload and participation mode on failure detection were 
separated. The participation mode was shown to have a dominant effect on the 
failure detection performance, with a failure in a monitored (coupled) axis 
being detected significantly faster than a comparable failure in a manually- 
controlled axis. 


INTRODUCTION 


In the last decade, a great deal of thought has been given to Category III 
landings and their implications. One area of intensive investigation centers 
around the role of the crew during the approach, and current thought is pola- 
rized around two extremes: 

a. The crew is in the control loop and flies the aircraft in accordance 
with instrument- generated steering signals. 

b. Steering signals are coupled directly into the autopilot, with the 
crew monitoring the system. 
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It is axiomatic that a pilot should be capable of detecting and identifying 
failures in the automatic landing system accurately, reliably and with minimal 
time delay. To this end, extensive studies have been conducted in which the 
pilot was treated as a controlling element in a one-dimensional task; his 
decision processes (Schrenk, 1969) and his adaptive behavior following a sud- 
den, change in the controlled plant dynamics were investigated (Young alt 
1964; Bhatak and Bekey, 1969). Other studies investigated the failure-detec- 
tion performance treating the operator as a pure monitor (Gai and Curry, 1975). 
In reality, however, the pilot is faced with multi-axes, not single-axis, 
tasks; although models for interference among multiple control tasks have been 
derived (Levlson, 1970) , the interrelationships among simultaneous control and 
monitoring tasks are not yet well understood (Levison, 1971). 

Young Qt al {op. ait.) found that in single-axis tracking tasks the human 
operator's performance as a failure detector was better when he was in the 
control loop; simulated Category III landing studies, on the other hand, have 
shown that the pilot's failure detection performance deteriorated when he was 
faced with manual control task, compared to the monitoring mode (Vreuls et alt 
196B). \^hen faced with split-axis tasks, pilots' monitoring and decision 
making were impaired (Monroe et at, 1968) and they sometimes completely over- 
looked the occurrence of a failure, presumably because of the increased work- 
load associated with split-axis tasks (Gainer et al^ 1967). 

It has been recognized that when the the role of the human changes from 
monitoring to that of an active controller corresponding changes take place in 
his workload level (Ekstrom, 1962; Wewerinke) . However, in pilot-performance 
studies to date these effects were completely confomded. It is the primary 
purpose of this investigation to separate these effects and to document pilot 
performance during a Category III landing as a function of the particular 
control mode at different workload levels. We wished to isolate and identify 
the effects on performance due to the variations in the control mode alone - 
and hence, variations in the operator's mode of behavior - apart from the 
effects on performance due to the variations in the workload level* 
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METHOD 


As stated* the purpose of this research was .the study of the pilot's short- 
term decisions regarding performance assessment and failure monitoring. We 
wished to investigate the relationship between the pilot's ability to detect 
failures, his degree of participation in the control task and his over-all 
workload level. Also, we wished our findings to be applicable to the general 
population of pilots who fly low-visibility approaches in commercial jet trans- 
port' aircraft. To this end, this research consisted of an experimental 
investigation which was carried out in a static ground simulator and which 
utilized fifteen airline pilots as subjects. 

The simulation capability included the ADAGE AGT/30 digital graphics com- 
puter and a fixed-base cockpit simulator. A mathematical model has been 
developed of a large transport aircraft in the landing approach flight enve- 
lope; the actual flight data of a DC-8 were used in the equations of motion, 
and the various parameters were later refined following a series of flight 
tests by a senior airline capcain with considerable Boeing 707/123 experience. 
Kon-linear phenomena such as ground-effect and stalls have also been included. 

An integrated- cue flight director system has been designed for this simu- 
lator, providing the capability to land the simulated aircraft manually in 
zero-zero conditions in a relatively satisfactory manner. Also, a two-axis 
autopilot has. been incorporated into the simulation which was capable of 
intefcepting and tracking the Instrument Landing System (ILS) , in either axis 
or in both axes, to touchdoxra. We also had the capability to add wind dist- 
urbances to the simulation to induce different workload levels. The wind gusts 
were modelled as filtered white noise with a cutoff frequency of tt/B rad/sec. 

The mathematical model was programmed into the AGT/30 computer which was 
linked via multiplexer channels and sense lines to the cockpit simulator. The 
cockpit was a mock-up of the captain's crew station in a Boeing transport 
aircraft (Fig. 1). The windows were frosted to eliininate external visual ref- 
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The controls Included an operational, spring-centered control colunm with a 
control wheel and rudder pedals, as well as four throttles, flaps, speed-brake 
and landing-gear levers and flight-director and autopilot controls. 

Apart from engine instruments and marker-beacon lights the simulator was 
equipped with three CRT screens, mounted one each on the main instrument panel 
at the captain^s and the first officer's stations and one in place of the 
weather radar screen. The screens were driven simultaneously by the ADAOE 
computer and presented the six standard flight instruments Cfig* 2): Airspeed, 
attitude- fiight director indicator, altimeter, instantaneous vertical speed 
indicator, horizontal situation (HSI) and radio-magnetic (RMI) indicators, as 
well as a DME digital readout and glideslope deviation and course deviation 
needles. The CRT screens were driven by the computer at a rate of 24 frames 
per second which was sufficient to produce flicker-free images. The informa- 
tion was updated at a rate of S/second. 

To measure the pilot's workload, a "warning light"-type subsidiary task was 
selected for the research. It consisted of two small red lights mounted above 
each other outside the subject's peripheral vision field, and a rocker thumb 
switch mounted on the left horn of the control yoke. 

The lights provided the stimuli. During the run the upper or lower light, 
with equal probability, was lit at a random time for a maximum of two seconds. 
A correct response by the subject consisted of turning the light off by a 
proper motion of the rocker thumb switch. The program recorded the number of 
times that the subject responded correctly to the warning light ("hits") and 
his response time (latency) for each response. Incorrect responses by the 
pilot, that is, not responding to an illuminated light or activating the 
switch the wrong way, were also counted and labeled as "misses". 

A workload index was computed from these data as follows: 

a. As each stimulus was presented for a maximum of 2 seconds, the total 
response-time ratio RTR for both "hits" and "misses" was computed by 

cumulative latency (?] T, ) 

i ^ 


.• . RTR = 


Total number of stimuli x 2 sec 


( 1 ) 






f 


b. A miss-rate MR was computed by 

^ _ Number o£ stimuli missed 
Total number of stimuli 


c. A workload index WLX was then extracted using the best least-squares 
fit weighing coefficients 


WLX 


0.78 RTR + Q.626 MR 
0.73 + 0.626 


X 100 percent 


(3) 


This measure of workload has been shown (Spyker et at, 1971) to be correlated 
with physiological predictors of workload with a correlation coefficient 
p = 0.6A6, significant at the ? < 0.005 level. 

d. Finally, we wished to eliminate differences between subjects which may 
have been caused by different subjects assigning different relative 
priorities to the primary tracking task and the subsidiary task. To 
this end the workload index of each subject was normalized, that is, a 
workload index of zero was assigned to the approach which resulted in 
the lowest workload measure for each subject and a workload index of 
100 was assigned to the approach with the highest workload measure for 
the subject. The normalized workload index on approach i of subject j 
was then computed by 


Normalized XJLX. . 

13 


WLX. . - 


min 

i 


{WLX. .} 


T - T 


X 100 percent 


(4) 


Experimental Design 


The experimental variables to be investigated in this study were the pilot's 
participation level in the piloting task, the workload induced by the control , 
dynamics and by external disturbances, and the pilot's failure detection 
performance. 

The experiment involved four levels of control participation.: 
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a. "Passive mo ni coring", with autopilot coupling in all -axes, including 
autothrottle. 

b« "Yaw manual", with autopilot coupling in the pitch axis and auto- 
throttle coupled. 

c. "Pitch manual", with autopilot coupling in the yaw axis only. 

d. "Fully manual". 

There were three levels of wind disturbance; 

a. No wind. 

b. A 45° tailwind of 5 knots, gusting to 15 knots. 

c. A 45° tailwind of 10 knots, gusting to 30 knots. 

Three failure conditions were used: 

a. No failure. 

b. Failure in the yaw axis. In this condition the autopilot, if coupled, 
or the flight director would steer the airplane away from the localizer 
course to intercept and track a course parallel to the nominal path but 
translated by a distance corresponding to one dot deviation (1.25°) at 
the point of failure occurrence. This resulted in a one-dot angular 
error about 100 seconds after the initiation of the failure. This type 
of failure was chosen, rather than a runaway failure, as it was quite 

• subtle and therefore it provided a good measure of the limits of the 

pilot failure detection capability. 

c. Failure in the pitch axis, which resulted in a one-dot deviation (0.35° 
of angular error) approximately 30 seconds after the occurrence of the . 
failure. 
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Failures were presented only between the altitudes of 1800 and 800 feet; each 
approach was terminated either at touchdown or when a positive rate of climb 
has been established following the initiation of a go-around by the subject. 
The selection of the failure altitude was randomized, as was the selection of 
the direction of the failure (left-right in a lateral failure mode, up-down in 
a pitch failure mode). Workload levels and failure detection performance were 
investigated in separate experiments, to avoid possible contamination of 
failure detection data by the presence of a concomitant subsidiary task; the 
"no failure" condition was incorporated In the design so that the subjects 
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would not anticipate a failure on. each and every approach. 


RESULTS AND DISCUSSION 


'It seems clear from Figures 4 and 5 that the side-task scores were sensi- 
tive to variations both in the disturbance level and the participation mode. 
Indeed, analysis of variance under the hypothesis that the effects of the 
disturbance and of the participation were additive revealed that the varia- 
tions in workload scores as a function of participation mode were significant 
at the P « 0.01 level and as a function of the severity of the disturbance 
- at the P < 0.03 level,' 

There was, however, no significant difference between workloads at the two 
low disturbance levels, namely, calm air and a quartering wind of five knots, 
gusting to fifteen knots, it was assumed, and it was verified by pilots’ 
comments, that Che components of the wind parallel and normal to the final 
approach path, 3.5 knots gusting to 10.6 knots, were not strong enough to 
induce workload significantly higher than that induced by piloting the sim- 
ulated aircraft in calm air. Consequently, these two disturbance levels were 
combined in the analysis 'and the data were treated as if there were only two 
distinct disturbance levels, "low" and "high". 

An additive model was used in the regression of workload scores on the 
disturbance levels and participation modes, to yield 

WLX(P,D) = W^CP) + W2(D) (5) 

where WLX is the normalized workload score 

P is' the participation mode 

D is the disturbance level 

18.7 for the fully -automatic mode 
36.6 for split-axis, yaw manual mode 
61.0 for split-axis, pitch manual mode 
72.9 for the fully manual mode 

, u ('n's ~ ^be "low" disturbance level 

an ‘ 2^ ~ 9.82 for the "high" disturbance level 


Ml(P) = 


NORMALIZED WORK LOAD 



Hormalized Workload Index at Four 

Participation Modes 

PI - Fully Automatic 

P2 - Split Axis, Yaw Manual 

P3 - Split Axis,. Pitch Manual 

P4 “ Fully Manual 


Figure A 


■ normalized work load 



Normalized Workload Index at Two 
Disturbance Levels 
Dl - Calm Air 

D2 r- 10 kt. Wind, Gusting to 3t> kt. 


Figure 5 


These values yielded workload-participation mode correlation significant at 
P < 0.001 and workload-disturbance correlation significant at P < 0.05. 

Detection performance was analyzed in terms of detection time and accuracy. 
Detection time was defined as the elapsed time between the occurrence of a 
failure and the verbal report by the subject that the failure has been detec- 
ted and identified. Accuracy was measured by the fraction of failures that 
were missed altogether. We differentiated between approaches in which a failure 
went unreported but which resulted in a successful touchdown and approaches in 
which a failure was missed and which did not terminate in a successful landing 
because of gross error in the failed axis. The latter are shown in Tables 1 
and 2; the numbers in parentheses represent the fraction of all missed 
failures, whether or not they resulted in a successful landing. 

In all, 90 approaches were flown in which a longitudinal failure occurred; 
of these, 8 went unreported, 6 of which did not terminate in a successful 
landing. Of the 90 lateral failures presented, 9 were missed; of these, 6 did 
not terminate in a successful landing. 

A very interesting pattern is obvious from Tables 1 and 2 and from Figures 
6 and 7: All failures in an automatically-controlled axis were detected in 
consistently short times; between 9 and 17 percent of the failures which 
occurred in a manually-controlled axis were not detected at all, and the ones 
that were required considerably longer detection times. The difference between 
the mean detection times in an automatic and manual mode was highly significant 
at Che P < 0.01 level. 

.We hypothesized that this difference in detection performance was due, in 
part, to the increased involvement of the pilot in the control task in the 
manual mode and, in part, to the increased workload levels associated with 
manual control; we set out to separate the individual effects of these factors 
on the failure detection performance. 

In Figures 6 and 7 the mean detection times of pitch and yaw failures, 
respectively, are plotted as functions of the corresponding mean workload 
levels for the four participation modes. The following relationships are 
evident: 

1. Detection times in a manually- controlled axis are longer than detection 
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TABLE 1 

Fraction of Missed Longitudinal Failures 
in Percent 


Participation 

Mode 

disturbance Level 

Overall 

1 

2 

3 

Monitor 

0. 

0. 

0. 

. 

0. 

Control Yaw 

0. 

0. 

0. 

0. 

Control Pitch 

12.5 

(12.5) 

0. 

(14.3) 

12.5- 

(12.5) 

8.7 

(13.0) 

Manual Control 

0. 

(12.5) 

14.3 

(14.3) 

37.5 

(37.5) 

17.4 

(21.7) 

Overall 

3.3 

(6.7) 

3.3 

(6.7) 

13.3 

(13.3) 

€.7 

(8,9) 


j 

i 


163 

















DETECTION TIME, SECONDS 


70 


O automatic 

^ MANUAL 
a SPLIT AXIS 


h 60 


30 


' 40 


5 i 


1 


■ i 


10 


20 


10 40 50 SO 

normalized work load 


70 80. 


90 


J'xgure 6, Mean Longitudinal Detection Times 


70 


h 60 


i 


so 


40 


Q AUTOMATIC 

A MANUAL 
□I SPLIT AXIS 


dk 

1 

i 


4 

2 


I 

1 


i 


Ik 


30 


- 1 * 


10 


20 


30 


40 50 60 

NORMALIZED WORK LOAD 


70 


BO 


go 


Figure 7. Mean Lateral Detection Times 


_J 

loo 


__I 

10 0 


165 


•times In an automatically-con trolled axis, 

2 , Detection times for lateral failures are significantly longer than 
detection times for longitudinal failures at comparable workload levels. 

3. Detection times increase in direct relationship to workload (p « 0.322 
for n=163 pairs). 

We assumed that the failure detection mechanism of the human operator acts 
similarly in both lateral and longitudinal axes; any difference in performance 
between these axes is due to differences in the plant dynamics and in display 
variables only, not to differences in processes internal to the operator- This 
assumption of equivalence between the lateral and longitudinal axes has been 
made, either explicitly or implicitly, by many investigators. It is based on 
the theory that the human behaves optimally with respect to his task (c£. 
Smallwood, 1967) in all axes, and that the operator adjusts his describing 
function to match the task (Young, 1969). 

Longitudinal and lateral failure- detection data were thus pooled; detection 
times were regressed on the type of failure (longitudinal or lateral) and on 
the control mode in the failed axis, with the workload index as a covariate, 
based on the following additive model: 


T 

detection 


A solution was 


Tq •f a(contol mode) + &(failed axis) + y(workload) 
obtained for the regression coefficients a, 3 and y; 


T 

detection 

where 1 

M “ 

0 


20.9 + 16.5 M +15.4 A + 0.10 WLX 

if the failed axis is controlled manually 
otherwise 


( 6 ) 

(7) 


and 


1 if the failure occurs in the 

A ^ 

0 if the failure occurs in the 
WLX “* the normalized workload index 


^ is measured in seconds, 

detection 


lateral axis 
longitudinal axis 


The relationship is plotted in Figures 8 and 9 for longitudinal and lateral 














failures, respectively. Mean, detection times at the corresponding mean work- 
load levels are also shoT.m for comparison. The model correlates well with the 
data, with p ” 0*531, significant at the P « 0.001 level. 


CONCLUDING REMARKS 


Our goal in this research was to identify the participation mode and work- 
load level which optimize the pilot’s failure detection performance; this 
subject is treated in considerably more detail elsewhere (Ephrath, 1975) • 

Our results indicate quite clearly that a coupled, fully-automatic landing 
with the lowest possible workload is called for in Category III operations, 
with the crew monitoring the progress of the approach via cockpit displays: 
Failure-detection performance in all other control modes was unacceptable for 
commercial operations. Performance monitors and fault annunciators may 
alleviate the problem somewhat but they are inadequate at altitudes below 100 
feet (Vreuls et al, 1968); also, they are not infallible, and additional 
warning lights and buzzers in the cockpit provide more opportunities for 
malfunctions and for crew confusion. 
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CHANGES IN PILOT WORKLOAD DURING AN INSTRUMENT LANDING 
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Introduct 

It 1 . j been, recognized that when the role of the pilot 
changes from that of an active controller to that of a monitor 
of an automatic system corresponding changes take place in his 
workload level (Ekstrom, 1962). However;, in pilot performance 
studies to date these effects were averaged over the entire 
range of a particular flight phase - usually the landing 
approach. Useful information may be lost by such averaging. 

Pilot's performance has been shown to vary greatly among 
different segments of an approach while flying by reference to 
instruments (DeCelles et alj 1970; Gainer et al, 1967; Ephrath, 
1975) . Such variations in performance may be caused by corres- 
ponding variations in the pilot's Instantaneous workload level 
during the landing approach. 

If indeed such is the case, and the observed deterioration 
in the pilot’s instrument- flying performance is caused by excessive 
workload levels, than it is evident that improved piloting per- 
formance may be achieved through cockpit designs that eliminate 
the causes of these excessive workloads. The purpose of our work 



is hence threefold: to examine the effect on workload of 

different levels of automation, to document the variations 
in the instantaneous workload level of a pilot during an 
instrument approach for a given level of automation; and to 
determine the particular phases of the approach during which 
the workload levels are highest as a prerequisite to identi- 
fying their causes. 

Method 


In the course of this study, fifteen professional airline 
pilots flew a total of 130 instrument landing approaches in a 
fixed-base cockpit simulator which was programmed to duplicate 
the dynamics of a large jet transport aircraft; the entire 
experiment was controlled on-line by an ADAGE AGT/30 digital 
computer . 

The pilots instantaneous workload level was measured by 
means of anon-loading warning-lights-type subsidiary task using 
two small Cl/8" diameter) red lights mounted above each other 
outside the pilot’s peripheral field of vision at 75° to the 
right of the center of the instrument panel. At random times, 
during the approach, the upper or lower light was lit with 
equal probability. A correct response form the subject con- 
sisted of a proper motion of a control-yoke-mounted, three- 
position switch; that is, pushing the switch up or down in 
response to the top or bottom light respectively. A correct 
response by the subject caused the light to turn, off; a "hit" 
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was counted and the subject's response time recorded. In the 
absence of a correct response the light remained on for two 
seconds, then was turned off, and a "miss” was counted. An 
incorrect response by the subject (that is, pushing the switch 
the wrong way} was also counted as a "miss”. 

After the light was turned off, a time delay followed, 
uniformly distributed between 0.5 and 5.0 seconds, and the 
process was repeated. The workload index, WLX, was computed by 

WLX = ion accumulative response time) -b S (number of "misses”) 

C2arf 3) (number of "hits" 5- number of "misses") 

which resulted In a workload-index value between 0 and 100. The 
coefficients a and 3 were chosen to maximize the correlation 
between this index and physiological measures of workload (Spyker 
et al, 1971) . 

A more detailed discussion of this subsidiary task and of 
the design of the simulator system can be found elsewhere (Ephrath, 
1975). 

Each subject flew twelve simulated approaches from a point 
well beyond the outer marker to touchdown at four levels of auto- 
mation or participation; fully manual, split-axis with manual 
yaw (pitch and power being controlled by an autopilot), split-axis 
with manual pitch (yax<r being controlled by an autopilot) , and 
fully automatic, with the pilot monitoring the instruments. Three 
approaches were flown by the subjects at each automation level, 
with the order of presentation randomized to eliminate learning 
effects. To this end, each pilot was also allowed ample practice 
in the simulator prior to commencement of the formal experiment. 
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Results 


Workload data were segregated by participation modes; 
mean scores were extracted at 300 feet altitude intervals 
between the altitudes of 2000 and 500 feet, and at 100 feet 
intervals between 500 feet and touchdown. The results, 
averaged over all pilots, are shown in Figures 1 through 4. 


Insert Figures Here 


Figures 1 through 4 reveal the asymptotic behavior of 
the pilot's workload as Inferred from the index. The work- 
load index is essentially constant at altitudes above 
approximately 500 feet, but there is a very marked increase 
in the workload at lower altitudes where the workload index 
appears to be inversely related to distance-to-go . This 
increase may be at least partially due to the non-linear 
increase in the display sensitivity. As the standard flight 
instruments display angular rather than linear deviations 
from the localizer course and from the glide path, the display 
increases in sensitivity as the distance to touchdown de- 
creases. However, this increase' in pilot workload was 
observed even when the simulated aircraft was controlled 
by the autopilot and the subject acted merely as a monitoring 
element; under autopilot control, the indicator instruments 
barely deviated from the zero-error position and, therefore, 
any changes in their sensitivity could not affect the monitor— 
ing subject's workload. This seems to suggest the effects 


of other factors such as, possibly, changes in the pilots' 
mental state arising from their awareness of the proximity 
of the ground and of the associated increase in the penalty 
for error. 

The observed increase in pilot's workload at low 
altitudes is quite interesting, as it is precisely this phase 
of a landing during which the pilot's instrument-flying 
performance becomes unacceptable and a transition to visual 
flight is mandatory; if visual reference with the ground 
cannot he established, the approach must be aborted (BeCelles 
et al, 1970). 

It is conceivable that display design might reduce the 
additional workloads imposed on the pilot as a result of 
increased display sensitivity; further studies are necessary 
however, to determine whether other factors are inherent in 
this situation and whether an improvement is possible in the 
pilot's workload - and his flying performance — during the 
last phase of the landing approach. 
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ABSTRACT 

A model of the pilot as a monitor of instrument failures during automatic 
landing is proposed. The failure detection model of the pilot consists of two 
stages: a linear estimator (Kalman Filter) and a decision mechanism which is 
based on sequential analysis. The filter equations are derived from a simpli- 
fied version of the linearized dynamics of the airplane and the control loop. 
The perceptual observation noise is modelled to include the effects of the 
partition of attention among the several instruments. The final result is a 
simple model consisting of a high pass filter to produce the observation resi- 
duals, and a decision function which is a pure integration of the residuals 
minus a bias term 

The dynamics of a Boeing 707 were used to simulate the fully coupled final 
approach in a fixed base simulator which also included failures in the airspeed , 
glideslope, and localizer indicators. Subjects monitored the approaches and 
detected the failures; their performance was compared with Che predictions 
of the model with good agreement between the experimental data and the model. 


INTRODUCTION 

The introduction of the "all weather" automatic landing system changes 
the role of the pilot during landing. Under normal conditions, the pilot is 
not in Che control loop, but his main task is to monitor the proper operation 
of the automatic system. This, of course, shifts his role from manual con- 
troller to decision maker. 

The problem of modelling the pilot as a controller has been addressed by 
several researchers, and satisfactory models exist using classical control 
theory (1) or optimal control theory (2). Models for the pilot as a failure 
detector have only recently been addressed (3), and some conjecture has been 
suggested (4). 

In this paper, a model of a pilot as a monitor of system failures is pro- 
posed and applied to an automatic landing. The model consists of two stages; 
a linear estimator and a decision mechanism. The linear estimator is the Kal- 
man Filter which is similar to that used in the optimal controller model; the 
outputs used here are the measurement residuals rather than the estimates. The 
decision mechanism is based on sequential analysis (5) modified for the spec- 
ial case of failure detection (6). Experiments were designea to validate the 
proposed model in which subjects monitored failures In simulated automatic ILS 
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approaches In a fixed-based jet transport cockpit. The results of these exp 
eriments are then compared to the prediction of the model. 


PROBLEM STATEMENT AND SIMPLIFICATION 


A functional block diagram of the failure detection model is shown in 
Figure 1. A basic assumption in the structure of the first stage (the estima- 
tion) is that the dynamical characteristics of the system that produces the 
input signals are known by the observer. Therefore, before the modelling of 
the failure detection system, we will discuss the model that the pilot has In 
mind (the internal model) for the airplane dynamics and control loops. 

The true airplane dynamics, when angular accelerations are neglected, can 
be defined by nine first order nonlinear differential equations. Two decoupled 
autopilots are used to regulate the vertical error between the aircraft posi- 
tion and the glideslope beam, and the horizontal (angular) error between air- 
craft position and localizer beam. In addition, a third control loop controls 
the aircraft airspeed. This configuration was used in the simulation that 
automatically landed the Boeing 707 dynamics used in our experiments (7). 

Since the pilot is outside the control loop his inputs consist only of 
the displayed variables on his instrument panel. If the control system is de- 
signed properly, these displays will show nominal values with variation due to 
outside perturbations. It seems reasonable to assume that the pilot will use 
the linearized version of the automatic system around the nominal values. In 
addition, we will assume that the longitudinal and lateral dynamics are de- 
coupled in the pilot’s model. The block diagrams of the three control loops 
are shown in Figure 2, and the basic configuration was taken from reference 8. 
Therefore the three closed loop transfer functions are given by 


5u lOCS-fQ.l) 

6u “ (S+8.a)(S+.98)(S4-0.13) 

^ where: 

6y 31.3 

66 (S+0. 5) (’s+5 . 5)’(S^+5 .’45S+11, 4) 

n 

5ij? _ 47 

6li/ “ l:s^+11^5S) (5^+1.53+0.81) 


( 1 ) 

u - velocity 
y - flight path angle 
0 - pitch (2) 

ijj - heading 

(3) 


The letter 6 is used to identify the inputs as perturbations rather than com- 
mands, and the outputs are the responses to these perturbations , The sub- 
script n is used because the input perturbations are modelled as zero mean 
white Gaussian processes. The source for these perturbations is usually the 
wind gusts, and therefore the inputs to the subsystem are correlated and are 
derived from the amplitude and direction of these gusts. Two of the above 
subsystems are of fourth order and one is of third order. Another integration 
of each subsystem output is needed to obtain the aircraft position. There- 
fore, we assume that the pilot uses only the dominant poles, namely the ones 


with the longer time constants. The final model that is used in the pilot 
model is 

- l/CS+l) (4) 

n 

» l/((S+0.5)(S+2.7)) <5) 

n 

U « 1/Cs2+1.5S+0.81) ( 6 ) 

'^n 

Note that the steady state gain and the steady state variance of the response 
to a stationary random input were kept the same as in the original systems. 

Having three decoupled systems, we can define 8 state variables by trans- 
forming equations (4) to ( 6 ) to their state space format. Define 
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The perturbations in the aircraft position in terms of the above state variables 
are given by (using the fact that Yo is small): 

6x = cos^fjQXi - vosinij^oxe 

6y = slntpoxi - vocosif'oxe (9) 

6 z = YoXi + V 0 X 3 

where the subscript 0 designates nominal values. 

All the variables that are displayed to the pilot can now be represented 
as linear functions of the state variables ^ 

1. Glideslope Indicator: yi = (-cosiiJo/x^+Yo/Xjj)xi+voX 3 /Xj^+vosiai|JoX 6 /x^ 

2. Localizer: yz =• Icosifio/ (1. 23-Xjj) ^+sln4Jo/ (!• 23-Xj^) ]xi 

+ [vocosij*o/Cl*23-Xjj)-vosini|io/(l»23-Xjj)^lx6 

3. Airspeed indicator: ye " xa 


I 


(pitch) 

(roll angle) 




4. Altitude indicator: 74 = xi, + xs/0.585 

yi* ■= VQXe/g 

5. Horizontal situation display: y$ » x? 

6 . Altimeter; y? =* Yo xi+ V0X3 

7. Vertical speed indicator: yg =ycSa+ vqXj, 

in the above equations is the nominal value In the x direction which Is 
me varying. 


XHE FAILURE DETECTION MODEL 

In the last section we suggested a simplified linear model which the pilot 
is assumed to use as a model for the instrument output dynamics. On the basis 
of these dynamics the failure detection model that was suggested by Gai (6) is 
used. As mentioned before, the model consists of two stages: a linear esti- 
mator and a decision mechanism- The linear estimator, the Kalman Filter, is 
shown in Figure 3. It should be noted that due to the time dependency of the 
measurements the Kalman gain K(t) is time varying. The filter produces esti- 
mates for the state x(t) and the measurements y(t) as well as the measurement 
error (residual) e(t) which is defined as 

£(t) = y(L) - yCt) (10) 

Any of these three quantities can be used as an Input to the decision mechanism. 
The observation residual is preferred for the following reasons; 

1. The state variables are non-unique variables that can be defined in 
many ways, while the observation residual is unique and well-defined for 
the subject. 

2. The dimension of the state is in general larger than the dimension of 
the residual. 

3. The residual is more sensitive than the observation to the effect of 
the failure (9). 

4. The residual is a zero mean white process (10) in the unfailed mode 
so successive observations are independent for Gaussian processes. 

In the instrument failure mode, we will assume chat a deterministic mean is 
added to the measurement so that the residual is still white Gaussian with the 
same variance but with a non-zero mean. 

Since the pilot is using 3 instrviments the problem of sharing of attention 
has to be accounted for. This Is done through the measurement noise in the 
observer model (11). If the pilot is observing more than one instrument, his 
observation noise for each observation is increased by a constant factor that 
is inversely proportional to the fraction of attention that he spends monitor- 
ing that specific instrument. Finally, it should be noted that although the 
state equations are decoupled, the Kalman Filter is a coupled 8 dimensional 
system because of the coupling through the measurements. The model of the 
estimation scheme Is shown in Figure 4. 
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The decision raechanism is based on sequential analysis (5). The classical 
sequential analysis uses the likelihood ratio i(m) as a decision function after 
m observations. Two criteria levels, A and B, are chosen, and the decision 
rule is 

if ACm) >_ A choose "failure” 

if Z(m) ^ B choose "normal" 

if B < Jl(ni) ^ A take another observation 

A and B are determined by the desired probability of false alarm P(FA) and the 
probability of miss P(MS) follows (5) 


A = (1-P(MS))/P(FA) 


B = P(MS)/(1-P(FA)) 


Since in our case, both distributions are white Gaussian with equal variances 
and means zero and 0i (failure), the decision function (for Gj >0) is (6) 

ra 

X(m) = I U. - 0i/2} (12] 

i=l ^ 


the upper criterion level is 


(lnA)/0; 


and Che lower criterion level is (lnB)/0i (14) 

The classical theory camot be applied directly to the failure detection 
problem because the basic assumption in the derivation was that the same mode 
exists during the entire period. A failure detection problem is characterized 
by a transition from the normal mode to the failure mode at some random time 
t-. In order to overcome this difficulty, we followed Chain (12) by: 

1. Resetting the decision function to zero whenever 3i!(m) is negative. 

2. Using only an upper criterion level Ai which is modified to keep the 
same mean time between two false alarms as before including the 
resetting. 

The value of Aj is related to A and B in equation (11) by the equation 
Ai - InAi - 1 = -{InA + (A ^l)lnB/ (1-B) } 

The modified decision function is shown in Figure 5 and the block diagram 
of the basic decision mechanism is shown in Figure 6. For the case 0i<O, the 
decision function is 

m 

X(m) = Z (e 4- 6i/2) (15) 

i=l 

and only the lower criterion level xs used. This criterion level Is 

-(InA )/0i (16) 

The final block diagram of the decision mechanism is shown in Figure 7. 

The operation of the proposed model is actually quite simple in principle. 
Its basic properties are; 

1. A high pass filter as a first stage to obtain the residuals 

2, Integration of the residual and comparison to a fixed threshold as a 
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decision mechanism (6). 

3. Only three parameters control the performance of the model: 

a. The parameter designating the mean of a "failed" process, 6i 

b. The signal to noise ratio of the observation noise in the 
Kalman Filter 

c. The probabilities of the two types of error P(FA) and P(MS), 

EXPERIMENTAL VALIDATION 
Method 

The Adage Model 30 Graphics Computer was used to simulate the non-linear 
dynamics of a Boeing 707 and the control loops (7). The output variables were 
fed into the instrument panel of a fixed based Boeing 707 simulator. 

The simulation included the last five minutes of flight prior to touch 
down from 10 miles from threshold and at 2500 feet altitude; the approach and 
landing were fully automatic. The failures that were used were instrument 
failures, so that they affected only the output variables and were not fed back 
CO Che system. In order to limit the experimental requirements, we considered 
only failures in two instruments, the glideslope (CS) indicator and the airspeed 
(AS) indicator. Four levels of failures were included for each of the two 
instruments. All failures were deterministic step changes that were fed into 
the instrament through a single pole low pass filter with 0.1 second time con- 
stant. The magnitude of the failures for the AS indicator were 

Cl = 2o ca >= 3a ca = 4a, = 5a (17) 

^ V V V v 

and for the GS indicator 

c. - Ogj c, . 1.50j,g c, = 2dg5 c. = 2.50j,j (18) 

where a and a ^ are Che standard deviations of the perturbation from Che nom- 
inal of^'^the displayed variable on the AS and GS indicators respectively. The 
cutoff frequency of these perturbations were ir/b radians/sec. Two random num- 
ber generators were used to choose the failure in each run; one determined 
the instrument and the other the size of the failure. In addition, a third 
random number generator was used to determine the time of failure, 

There was a single failure in 90% of the runs. The high percentage of 
runs with failures was chosen to provide enough data in a reasonable experi- 
mental time. There was no feedback, to the pilot concerning his performance 
because it was felt that such feedback, would bias his decision, by driving him 
to try to compensate for previous errors, or to overrelax after several correct 
decisions. 

Each subject participated in three experimental sessions each of which 
included 16 runs, or a total of 48 runs. When the pilot detected a failure, 
he pressed a button and the run was terminated. Otherwise, the run would go 
until touch dovm. After termination of each run the subject was asked to fill 
out a form in which he stated which instrument failed and how he detected the 
failure. 
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At the beginning of each session, a set of Instructions was read to the 
subject. In particular, he was told that failure would either be in the AS 
or GS Indicator, but that he could use other instruments for verification. 

Results 

The experimental results for two subjects are summarized in Table 1. The 
tabel shows the mean and standard deviation of the detection time for failures 
in the two instruments. The results are also shown in Figures 8 through 11. 
These figures include the mean detection times that were predicted by the model. 
These mean values were obtained by the use of a Monte Carlo simulation, with 
the same number of runs as in the experiment. The values for the three model 
parameters were: 

SNR = 36 P(FA) == P(MS) = 0.05 = 0.25 (19) 

The level of P(FA) was determined on the basis of the actual false alarm rate 
that was found in the experimental data. For both subjects, the predicted 
results seem to fit the experimental data well. Of course, better fit could 
be obtained by change of the parameters in equation (19). 


CONCLUSIONS 

In this paper we proposed a model for the performance of a pilot as a 
failure detector of instrument failures during an automatic landing. The 
model consists of two stages: a linear estimator and a decision mechanism. 

The linear estimator is the Kalman Filter which is determined from a simplified 
model of displayed-variable dynamics that are used by the pilot. The filter 
also accounts for the pilot's time sharing between instruments through the 
observation noise. The decision mechanism is based on classical sequential 
analysis with some modification for the failure detection case. 

An experiment designed to test the validity of the model is described. In 
this experiment, subjects had to detect failures in the glideslope and air- 
speed indicators during a simulated landing in a Boeiiig 707 fixed base simu- 
lator. The results show that the predicted detection times fit the experi- 
mental data well. It should be remembered, however, that only instrument 
failures in the form of a change in the mean were discussed. Additional 
work is needed to verify the model to Include failures that involve dynamic 
changes as well. 
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FIGURE 4. FUNCTIONAL DIAGRAM OF FAILURE DETECTION DURING AUTOMATIC LANDING 
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INTRODUCTION 

Despite a vast accumulation of operational experience with the conduct of 
low visibility instrument approaches, little is understood about the decision- 
making behavior of pilots who fly these approaches. Likewise, there is little 
information regarding the man, system, and task-related factors which influence 
this decision-making behavior. Such information is essential for the rational 
design of new systems, or for the redesign of existing systems in order to 
correct knoiim deficiencies. 

A major problem which has Inhibited the study of pilot decision-making be- 
havior in the laboratory has been the unavailability of tasks which incorporate: 
the essential cognitive features of the real task, and which include those 
motivational or stress-inducing features known to influence decision-making 
performance. This paper describes a task which was designed to simulate the 
cognitive features of low visibility instrument approaches, and to produce, 
controlled amounts of subjective stress in pilots serving as subjects in exp- 
eriments using the task. 

Pilot behavior during low visibility instrument approaches can be analyzed 
Into at least two major categories: one is the continuous closed- loop manual 

tracking behavior necessary to control the aircraft, and the other is the cog- 
nitive, decision-making behavior required to make the decision to continue the 
approach and landing, or to execute a missed -approach. It is the second cat- 
egory of behavior with which we are concerned here. 

The difficulty of a decision-making task is, in part, determined bv the 
uncertainty of the data used to make a decision. For example, the decision to 
"go around" is a relatively easy one if, at the missed approach point, there 
is nothing to be seen outside the aircraft, or if the approach light.s and 
runway have been clearly visible for the last two miles of the approach. It Is 
when the approach lights or runway are barely visible, and tiien only intermit- 
tently, that the decision-making task becomes more difficult, 

A second class of variables which are known to influence the outcome of 
decision-making tasks is best illustrated by the various kinds of osycliological. 
stressors acting upon the pilot. Of particular interest here are ihe pressures 
perceived by the pilot to complete the approach, to make an on-time arrival, 
to save fuel, and even to save "face". 

f 

We have assumed that it is necessary to use a simulation task which incor- 
porates both kinds of variables, informational and psychological, to success- 
fully study pilot decision-making behavior in the laboratory. The ta.sk below 
was designed to meet those requirements. This paper describes our preliminary 
experiments in the measurement of decisions and the inducement of stress in 
simulated low visibility approaches. 
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METHOD 


A schematic of the apparatus as seen by the pilot-subject is shown in Fig. 
1. The buttons available to the subject are RVR (to request an RVR reading), 
turn rate buttons (left, 0, or right), and GA, the go around button to initiate 
a missed approach. 


In the central portion of the CRT is a plan view of the approach. In the 
lower part of the screen are three dots corresponding to the position of the 
approaching aircraft (present position, position one second ago, and position 
5 seconds ago). In the center of the screen are two pairs of dots correspond- 
ing to the middle marker location, equivalent to the 200 foot decision height 
for a Category I approach. Farther up the screen are the runway outline, 
threshold, and three pairs of approach lights or lead-in lights. .-\bove that 
are scores posted for the results of any one trial: on this approach the sub- 

ject would receive 100 points for a safe landing and -40 points for a missed 
approach. On the left of the screen is an RVR scale with two indices corres- 
ponding to 0 RVR and that for the legal minimum (2400 feet). On the right 


side of the screen is an altimeter which has a dynamic range of 0 to 220 feet. 
The pointer indicating altitude is pegged at the upper right until tne ai .''craft 
nears the middle marker; as the aircraft passes through the middle marker, the 
indicated altitude passes through 200 feet- 


A random wind disturbance from the side (correlation time ol 50 seconds) is 
introduced to provide a moderately easy control task for the pilot. Control is 
maintained by pushing one of the three turn-rate buttons. The aircraft has 
the capability of being in either the 0 turn rate (constant heading) or a stan- 
dard turn rate to the right or left. The pilot's task in these approaches is 
to "fly" the aircraft through the gate, over the approach lights, and on to 
the runway. (The aircraft's position sho;vn in Figure lis close to the initial 
condition point.) Only lateral position is important, for if the pilot crosses 
the extended threshold line, but is not over the runway a crash is recorded. 

If at time before the aircraft crosses the extended thres'nold Jlne, cho 

pilot hits the go-around button, a standard rate left turn is initiated until 
the heading reaches 6C‘^ from "North" at which time the computer program assumes 
that a missed approach was made. 

The runway and approach lights may appear either to the riglit or Ic t L oi 
the middle marker center line, and may be closer or farther away than Clio nom- 
inal position to represent electronic guidance errors. Tliis is the approj'.r lace 
aircraft-centered view, and simulates the case when one is flying tliu I. i.S vicli 
needles exactly centered, but finds the runway to the left or rigiit when break- 
our occurs, and the case when one is either high or low of the indicated 
altitude. 

The slant range "visibility" is Included in cne program, even though the 
Intensity in the CRT has only two values (on, off). There are 5 "characters" 
drawn by the PDP-12 graphics system under visibility control: the tliree pairs 
of lead-in lights, and the right and left halves of the runway / thre.sho Id ji.-’.hcs. 
Should the center of any of these five characters be within a square (cciiC‘'red 
at the aircraft position) whose half-width is the slant range visibility, then 
this character will be turned "on" and will be visible. The approach lights 
are turned off as one gets close to each pair, to simulate their passing 
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underneath the nose of the airplane; this also prevents the subject from ob- 
taining unrealistic lateral guidance information. 

A computer program was written to generate files of approach trajectories 
and currently has a catalogue of nine approach trajectories. Five of these 
trajectories have constant (but different) slant range visibilities leading to 
the following effect: when the middle marker is passed, nothing is in view; 

soon the first approach light appears, followed by the second and then the third; 
as the first approach light is neared, it disappears (passes underneath), and 
theii the runway /threshold lights suddenly appear and a safe landing can be 
accomplished. The decreaoing slant range visibility in this group of five tra- 
jectories is such that one must proceed farther and farther beyond the middle 
marker (or below decision height) before the first approach light is sighted. 

The fifth of these five trajectories is zero-zero visibility, so the approach 
lights and runway /threshold lights never appear. The other four trajectories 
correspond to: 

(1) a high vi^ ility approach (runway s>nd approach lights are visible as 
shown in Figure 1 at all times) 

(2) an extremely optimistic RVR reading, but very low slant range 
visibility 

(3) passing through a fog bank after Initial acquisition of the approach 

lights: the approach lights and runway lights "drop out", only to 

reappear after three to four seconds 

(4) fog bank as in (3), but the approach and runway lights do not 
reappear, 

PROCEDURE 

In the first set of experiments, we had the following objectives: 

(1) to structure the experimental setting to make the pilot as aversive 
to a crash in the simulator as he would be in real life 

(2) to alter the decision strategies by manipulating the relative values 
of a landing and a missed approach. 

The first objective was desirable to make the decisions as meaningful as pos- 
sible. After "sacrificing" several pilots, we finally arrived at the following 
procedure , 

As the subject is led into the experimental chamber, he is shown a poster- 
sized list on the wall of people whc have previously been subjects in the exp- 
eriment. Each subject is listed by name, organization, and score (the total 
number of points accumulated over the 50 data trials). The first subject on 
the list was a fictitious one (in this case) and in place of his point score 
was the word CRASHED in bright red letters. The experimenter writes in the 
subject's name and organization (e. g. Joe Jones, TWA) and leaves the score 
column blank. The subject is told at that time that should he crash during 
the data trials, even if on the first data trial, his services are no longer 
required. That is, in terms of the experiment, he is "dead". 

It was obvious to the subject at this point that he was committed to follow 
through the experiment, and the idea that he might crash and have the event 
recorded for all to see had a very noticeable effect on almost all the subjects. 


Each of the pilots was allowed a total of 25 -practice approaches, the first 
10 of which were high visibility approaches so that he could become familiar 
with the dynamics of the simulation, the wind and turbulence levels, and the 
layout of the approa-.'h lights, runway, etc. After a brief rest period, the 
pilots participated in the 50 data trials — it was during these trials chat 
a crash would mean imniediate dismissal. 

The 50 data trials were composed of six "wild card" approaches, e.g. the 
approach lights dropping out and then reappearing, or approach lights dropping 
out and not coming back. The remaining A4 trials were the ones examined for 
pilot decision-making behavior and consisted of eleven replications of four 
meteorological visibility levels of 0, 20, 30 and 40 display units, where a 
visibility of 50 corresponded to having the first approach light come into 
view as decision height was reached. 

These 44 approaches were assigned go-around scores of 100 points for the 
highest visibility doira through -80 points for the lowest visibility and were 
not assigned randomly, but in a manner which we though would make the decision 
most difficult. In general, a negative score corresponded to a low visibility 
approach and a high positive score to the high visibility approaches. 

The data recorded during each approach consisted of a "frame" composed of 
the current x,y position, the displayed RVR, the state of the turn-rate control 
and the state of the go-around button. These frames were recorded whenever a 
control action was executed, an RVR request made, and when the go-around button 
was pressed. From these data, we can infer the number of times the RVR vv^as 
requested, the control activity, and the altitude and cross track error at the 
time of go-around should one be requested. 

QUESTIONNAIRE 


Thirteen pilot-subjects participated in the test and completed a question- 
naire, but as the simulation was changed after the first three pilots, they 
were not Included in the data regarding the simulation itself. Of the remaining 
10 subjects, 8 are airline pilots and 2 are IFR rated NASA employees. 

The questionnaire consisted of 3 major parts: recent experience in low 

visibility approaches and missed approaches; fidelity of Che decision simula- 
tion; and stress ratings for actual low visibility approaches and the simula- 
tion. 


Recent Experience 

Of the 10 pilots completing the questionnaire, 7 had made a total of 37 
Category I approaches within the last 12 months (six of these 37 approaches 
were military approaches). Only 2 missed approaches were made by these 7 
pilots. When asked what were the most common causes for executing a missed 
approach (based on their experience), the 3 most frequently mentioned items 
were 

runway alignment/crosswinds 7 times 

visibility 5 times 

other traffic 3 times 
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simulation Fidelity 


The subjects were asked to comment via the questionnaire about the simul- 
ator fidelity only with respect to the decision of whether or not to continue 
an approach . This was done both on a semantic differential scale (Totally 
Unlike - Completely Identical) and by soliciting comments on the similarities 
and dissimilarities of the simulation to an actual low visibility approach. 

The ratings of the subjects are shown in Table I where it is seen that the 
mean fidelity rating is 5.2 with a standard deviation of 1.87, indicating the 
usual dispersion in intersubject ratings. 

Comments on the similarities of the simulation to a low visibility approach 
detailed the assimilation of information through different sources (RVR, alti- 
tude, and runway alignment), l^hen commenting on the dissimilarities, 3 pilots 
mentioned the lack of danger ("one will not die if you miss", "... lacks the 
element of danger") . Two of the pilots mentioned that in a real approach more 
reliance would be placed on decision height, i.e., chat it is a cut and dried 
decision (a go, no-go situation). Another commented that he felt the reward 
structure was not correct because in actual flight the rewards for going 
below minima may be the loss of job, etc, whereas reward here is a higher 
point count. 

There were other comments made about dissimilarities of the simulator and 
the actual approach: three pilots mentioned that the visual cues were differ- 

ent, and one pilot mentioned the fixed turn rate characteristics of the simu- 
lator. Those were offered even though the question specifically ashed about 
the similarities of decision making; cither the questions were misunderstood 
or these factors really do influence the decision. In either case, we felt 
that these latter two factors are of secondary importance in the light of the 
other dissimilarities mentioned by the pilots. 

Stress Ratings 

The pilots were asked to rate the stress of the experimental tas '' and an 
actual low visibility approach on a semantic differential scale (N’ot at all 
stressful - Extremely stressful); the results are shown in the other columns 
of Table I. We have added columns showing the difference in stress rating, 
and the simulator stress (racing) as a fraction of the actual stress (rating). 
Of these 10 subjects, three felt that the simulator was at least as stressful 
as an actual low visibility approach. At the other extreme, is subject number 
six who reported that the simulator "lacks the element of danger". 

RESULTS 


One pilot misunderstood ’the instructions because he Initiated a missed 
approach after safely crossing the threshold many times during die 50 data 
trials and therefore received less than full point score. His data were not 
analyzed. 

Learning 

A statistical test was used to ascertain whether or not a learning effect 
was present for the subject group by performing an analysis of variance on the 
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pull-up altitude for those 11 approaches made under zero/zero visibility 
conditions. This analysis of variance included the approach number as a co- 
variate and, if significant, would suggest a linear trend in pull-up altitude 
with trial number. The results of this analysis of variance indicated that the 
covariate of approach (trial) number did not contribute a statistically signi- 
ficant linear component to the pull-up altitude. Thus the linear trend was 
ignored in the remainder of the analyses. 


Classification of Approaches into 
Land / Go- Around 

The 44 approaches for each of the 9 pilots were examined for classification 
into the classes of land /go- around using a stepwise discriminant analysis pro- 
gram (BMD 07M) . The variables included for the selection in the stepwise dis- 
crimination were the following: 




n 




n = 1,2,3 


1,2,3 


<V/Vj^)(S/S,^) 


Actual visibility on the approach 
(0 < V/V < 1) 

Score increment for selecting a go- 
around 
(-0.80 < 


3/S^il.O) 


Interaction between visibility 
score 


and 


LOG Visible localizer deviation at 

go-around 

The linear, quadratic, and cubic component of visibility and go-around score 
are self explanatory and the interaction term was included to test for its 
significance. The localizer deviation was also included and taken to be the 
maximum visible localizer deviation. It was set to 0 on the 0/0 approaches 
since it would not be available to the pilot, and was also set to zero on the 
approaches which were successfully completed. 

The significant variables selected by the stepwise discriminant program 
are shown in Table II for each subject; these coefficients have been normal- 
ized so that the coefficient of is unity. In this table, an incr...]se 

in the discriminant function will “ put the approacli into the "land" class 

The resulting classification using this discriminant function is also shown in 
the table and it gives very good results on these data (although on..- musL be 
cognizant that the classification is performed on the data from which the dis- 
criminanf- function was determined.) 

The major points to be ascertained from the table are first, the go-nround 
score was almost useless as a basis for discriminating among approaches with 
the exception of Subject number 6. (This particular discriminant funcLion 
must be treated with care since it incorporates almost ail the variables and 
includes a sign of V/V which is opposite from all the other discriminant 

functions.) The secontT point of interest is the nearly equal coef f i.c iemt on 
the quadratic component of visibility, indicating the decrease in effective- 
ness of actual visibility in classifying an approach into "land". Thus the 
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contribution of visibility to the discriminant function varies from 0 (for 

V/V„.„ of 0) to 0.4 (for V/V„.„ “ 1.0) for most subjects, 
nAJv nAA 

The coefficient for the localizer deviation can be used to determine the 
sensitivity of the cross-track error in classifying approaches into "land" or 
"go-around". This result is shown in the right hand column of Table II and 
Is the localizer error (in degrees) which has the same effect on the discrim- 
inant function as a full-range change in slant range visibility. This gives 
the Importance of localizer error relative to visibility in determining the 
classi f ication. 


A DYNAMIC DECISION MODEL 


In this section, we briefly describe a decision-theoretic approach to the 
modelling of pilot decisions during the simulation of low visibility approaches . 
A straightforward application of the Subjective Expected Utility (SEU) theory 
is complicated by the dynamic character of the decisions since the theory re- 
lates to static decision alternatives, rather than the everchanging situations 
experienced by pilots. Nonetheless, we have developed an extension (based 
SEU models) which appears to be plausible, and it is one which we think will 
be a valuable tool in further investigations. 

The dynamic decision model is based on the assumption of the existence of 
a decision function which can be written 

D(V,S,L,h) - (1) 

where D is the decision function, an explicit function of visibility (V), go- 
around score/incentive (S), the cross-track or localizer deviation (L) , and 
the current altitude (h) . This decision function is the comparison of the SEU 
for snd the SEU for making a missed approach or a go-around Under 

the assumption that Che utilities for landing, crashing, and going around are 
independent of the probabilities, Subjective Expected Utilities take the form 

^GA ^ AROUND) = U(S) (3) 

where P, AMn subjective probabilities for landing and crash- 

ing and^ U(*) IS the corresponding utility. These expressions show the 
dependence on the approach variables V, L, h and the incentive for going 
around S . 

A schematic plot of the decision function and how it might change with 
altitude is shown in Figure 2. We have displayed possible variations of 
during an approach and its comparison with 0^^ which remains constant through- 
out the approach. If at any time the SEU of ’ landing becotitcs less than that 
of going around, the decision is made to initiate a missed approach. 

A missed approach, denoted by the solid line of Figure 2, is a sketch of 
how the SEU of landing might behave during an approach for which the approach 
lights are never sighted. The SEU of landing decreases with altitude because 
the' subjective probability of landing is decreasing (and that of crashing is 
increasing) until, at point A, a missed approach is initiated. The SEU of an 
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approach shown by the dotted line starts out higher than the previous approach 
(perhaps because of a larger reported RVR) and it, too, decreases until point 
B when the approach lights are sighted. This causes an immediate jump in the 
probability of landing (hence the step change in SEU of landing) and from this 
point gradually increases until point C when the landing is successfully accom- 
plished. If, on the other hand, the aircraft starts deviating from the local- 
izer (say) at point D, and Che pilot has difficulty in stabilizing the approach, 
then the subjective probability for landing will decrease causing a decrease 
In the SEU for landing until point E where a missed approach is initiatied 
because of misalignment with the runway. 

The point of greatest interest, of course, is at the instant of deciding 
to initiate a missed approach. At this instant the SEU of landing and go- 
around are equal, i.e. 


where h* is the altitude at which the go-around decision was made and is 
written 


h* == g(V,L,S) (5) 

To test for the possible existance of such a relation, we performed a stepwise 
regression (using BMD 02R) on the go-around score » score/vis lb iiity 

Interaction localizer deviation (L) . The main 

effects of meteoroTogical visibility were not included because the majority of 
the go-around decisions were made under 0/0 visibility conditions. Table III 
shows the regression coefficients selected by the stepwise regression program; 
those coefficients vjhich have a non-zero value as indicated by Student's t 
test are indicated with an asterisk. Note that the data for subjects 2,3,4 
are not included because the stepwise regression program did not find a signi- 
ficant regression on the variables indicated. The multiple regression coeffi- 
cient is highest for those cases for which few go-around decisions were made 
during approaches when any visibility existed (recall that 11 of the approaches 
shown were made under 0/0 conditions). While the coefficients indicate the 
type of behavior one would expect, e.g. an increased decision altitude due to 
increased go-around scores, these data must be considered preliminary because 
of the experimental design (see the discussion section) . 

These results are quite encouraging and indicate that the subjectively 
expected utility model proposed here may lead to a valuable viewpoint from 
which to examine pilot decision making during low visibility approach. 

SUMMARY 

Stress: 

One of the major goals outlined for this preliminary set of experiments 
was to investigate methods of applying psychological stress analogous to the 
stress of an actual low visibility approach. It was found that the stress 
rating in the simulation, as reported on semantic differential scales, was an 
average of 0.8 times the stress rating of an actual low visibility approach. 
The success in applying the stress was not uniform, however, for several 
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subjects reported they "would not die" if they crashed in the simulation. 

Other subjects remarked that descending below the decision height of 200 feet 
would result in censure by regulatory or company authorities; this seemed as 
important to them as the prospect of an (unlikely) crash. Thus it may be 
meaningful to include another penalty in the simulation; if the subjects are 
"caught" descending below minimums, they will be penalized (say) by a score 
equivalent to two or three landings. 

Go-Around Incentives 

The range of go-around scores did not induce as much behavioral change as 
was expected, and the results of the experiment indicate that considerably moi-e 
score differential will be required to induce pilots to initiate a missed ap- 
proach. For example, when offered 100 points for the no-tisk go-around or 100 
points for a successful landing, most pilots initiated the approach and almost 
all. continued until touchdovm, even though a riskless go-around was available. 
This suggests that the utility of landing is strongly affected by the accom- 
plishment of this feat, or that the level of risk taking in a go-around is too 
low and this alternative does not present enough of a challenge to the pilots. 


Behavioral Models 

One model of behavior which may apply here is the Theory of Achievement 
Motivation (Atkinson, 1964), This model is of a form similar to the SEU model 
described above except that the utilities depend on the subjective prob.ibili- 
ties; for example, the utility of succeeding at an easy task is low, while 
succeeding on a difficult (iiigh risk) task is high. Conversely, the utility 
of failing on an easy task is low (one loses face), whereas there is no dis- 
utility (loss of face) on failing to succeed on a risky task. Risk taking 
behavior is said to be determined by two personality traits; need for achieve- 
ment and test anxiety. The Theory of Achievement Motivation predicts that 
those individuals with a high need for achievement and a low test anxiety will 
take an intermediate level of risk, whereas individuals with a low need for 
achievement and high test anxiety will take extreme levels of risk: a lov/ 
level of risk to insure success, or a high level of risk in which success is 
not really expected. Atkinson draws analogies to aspirations In employment as 
well as more quantitative behavioral tests such as tlie "ring toss" .axporimenC . 
We have conducted some informal experiments at MIT using a ball-toss iiaradigm 
involving two levels of difficulty; the results of this undergraduate student 
project will be reported elsewhere. 

A preliminary attempt was made to apply the Theory of Achievement Motiva- 
tion to the experimental results described above. An obvious measure of test 
anxiety is the stress ratio recorded by the subjects (stress in the s imuiation.'^ 
stress in actual low visibility approach) , Measures of success and need for 
achievement are ambiguous, hov;ever, since point score and number of landings 
may be considered as measures of both. 

Nonetheless, if one considers (a) the final score as a measure of success; 
(b) the stress fraction as the measure of test anxiety; and (c) the number o*^ 
landings as the measure of the need for achievement (e.g. sticking witli an 
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approach through approach light dropout, etc.)* then classification of the 
pilots on this basis (b) and (c) , indicates the following: the pilots with the 

lover success level (lower score) exhibited a low need for achievement and high 
test anxiety (as the theory predicts), but the pilots with a higher success 
level (higher score) exhrlhited not only a high need for achievement but a 
higher stress level (rather than the theoretically predicted low stress). 
.'Although there are not enough subjects to validate this conjecture statistic- 
ally, it may well be that the Theory of Achievement Motivation is not applicable 
in those cases where the result of the failure is catastrophic, and that mod- 
ification to the theory may be required for situations such as are considered 
here. 

In summary, both a dynamic version of Subjective Expected Utility Models 
and (a modified) Theory of Achievement Motivation may be useful in describing 
decision behavior of pilots in a simulated low visibilicy approach. 
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Figure 1 Subject's display panel 
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FIGURE 2. 


Decision fTinction versus altitude 


Subject/ 

Organisa-' 

tion 

Simulator 

Fidelity 

Rating 

Stress 

Rating 

-Actual 

Approach 

Stess 

Rating 

Simulator 

®sim"®act 

®SIM 

^ACT 

1/A 

7 

3 

2 

-1 

.67 

2/B 

7 

4 

6 

2 

1.50 

3/C 

3 

7 

3 

-4 

.43 

4/B 

3 

8 

5 

-3 

.62 

5/B 

7 

8 

6* 

-2 

.75 

6/A 

4 

7 

2 

-5 

.28 

7/A 

7 

6 

6 

0 

1.00 

8/C 

6.5 

5.5 

5.5 

0 

1.00 

9/D 

2.8 

7 

4 

-3 

.57 

10/D 

5 

8 

6,5 

“1.5 

.81 

Mean 

5.23 

6.35 

4.60 

-1.75 

.763 

S.D. 

1.87 

1.73 

1.73 

2.01 

.344 


* Indicated a change to 2 later in the trials 


TABLE I Semantic differential ratings of simulator fidelity 

and stress 
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Table Coefficients of Discriminant Function (Normalized such that 1 






























TABLE III Regression coefficients for altitude of go-around decision 
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Program Abstracts/Algorithms 


A multinomtal maximum 
likelihood program (MUNOML) 


RENWJCKE. CURRY 

Man- Vehicle Laboratory 
Department of A emnautics and Astronautics 
Massachusetts institute of Technology 
Cambridge, Massachusetts 02139 

In our research on moJeJing sensory and decision phenomena, 
we Were soon confronted with the task of evaluating both old 
and new models using both old and new data. Rather than design 
an ad-hoc estimation program for e.ich new niodet, as is typically 
done, we developed an "executive'’ program that provides a 
general method for estimating parameters and simultatteously 
provides flexibility for accommodating new models with a 
miniimmi amount of programming. Our experience with canned 
computer programs lias been cc|Uivocai, so we decided to provide 
only tfic general framework and let the user accomplish the 
objective of estimating parameters for his particular model by 
writing a new subroutine within the constraints of the executive 
program. 

The most comnion class of distributions for svhicli parameters 
must be extracted arit multinomial distributions resulting from a 
siirtioliis-response dassification, e.g.. binary responses (YbS-XO 
or two alternative forced chuiee methods), the method of 
successive categories (rating scales), or transition probabilities in 
a Markov chain. Aithougii a number of rnelliods exist !\ir 
osliniating sucti parameters (Restle, 1971) we have chosen the 
maximum likelihood method and have implemented the scoring 
of Rao to adjust the parameters from one iteration to the next. 
We have chosen the maximum likelihood (,\fL) method because 

(1) it is a member of the class of consistent asy mptottcaiiy 
normal estimators (CAN) under suitable condiiiorrs (Rao, 1973): 

(2) it will easily handle situations in whieii all I’le responses fall 
into one category; (3)t!ierc are many situations in which the 
maxinium likelihood estimator can be shown to yield iinicjue 
estimates for parameters (ocher estimation Icciiniciucs may liavo 
this property as well); and <4) it is tlie only one exhibiting first 
order efficiency (R.io, 1 973). 

Other authors have successfully used the inaxirmitii likelihood 
method wiili R.jo’s scoring mclliod in signal detection type 
models (c.g.. Uorfrnan and Alf. 196H'. Dorfman and Alt', 1969; 
Dorfnian, Deavers, it Saslinv, 1973; Grey & Morgan, 1972; 
Ogilvie & Crceliiian. 1968). I ids jirogram differs pnmanfy in its 
ability to estimate parameters (or a wide variety of niodel.s. 

A BRIRF Rt VII.W OF 
MAXIMUM LIKULIIIOOD liSTLMATJON 


j(j = 1, , . .jt) given stimulus i(i - I m). Thus, wc have a 

series of multinomial distributions (multimultinomial), and it is 
assumed that the responses in each of these distributions are 
independent. Under tfiesc conditions, the likelihood function for 
the multinomial distribution is proportional to L(x) (Rao, 
1973). 


tn n , 

L(x)= ri n P^U(x), (1) 

1=1 1=1 « 

where x is the S-diniensional parameter vector to be estimated 
and Ty is the number of "j” responses observed when the 
stimulus is specified at condition "i," 

As is usual, we will consider the rsatural logarithm of tlie 
likelihood function, since it is easier to deal witli numerically 
and analytically, and its maximum occurs at tJie same point as 
tJiat of the likelihood function. 

m ft 

En Lfx) = E(X) = .S ryCnPy(x). (2) 

The first order necessary condiiions for an ex'rcmurh poini of 
the likelihood function are given by the (components of the) 
vector eiiuation; 


1 


ap 




ax 


5 = 0 . 


(3) 


ax i=i j= 

ViTiether any extremum pornt s-itis lying (3) is a logal maximum, 
minimum, or saddle point is foiin.i by examining the eiacnvaiiics 
of the matrix of .tecond partial derivatives or the Hessian inatri.x 


B _ _v ^ ^ii iiPjj OPjy 1 

“ i=i j=ipg“'(x) ax (ax 


in ft *'ii S’Pjj 
^ i=i j?i 


(4) 


Finding the Maximum 

Sina' the necessary conditions for an extremum (3) arc rarely, 
if ever, analytically tractable, an iterative approach U called for. 
The general form of the parameter updating techniijiic usually 
takes the form 


Vi-Kl = -'i + . (5) 

3x 

wl'icre Is the parameter estimate at the iih iteration, av/3\ is 
ihe gradient evaluated at (hc' point Xj, anil the matrix Aj 
determines the algorithm by which one socks the niaximitm ol 
tile likeliliood function. One approacJi, Rao's scoring method, is 
to construct the matrix A from an approximation to the 
information matrix 


Necessary Ginditions for a Ma.ximum 

In this section, svo surnmari/e the pertinent aspects of the 
maximum likelihood (ML) estimators pertaining to tlie practical 
considerations of carrying out the estimation process. ISec tlie 
excellent descriptions by Rao tl973) and Kendall & Stuart 
(1967) for a more detailed discussjoii of ML estimators.] Wc will 
formulate the problem as it is usually encountered in behavioral 
research, that is. in the form of a multinomial distribution 
containing parameters to be estimatcij. Wc specify by Py the 
theoretical conditional proliahiliiy of observed response 

Sponsored by NASA Grant NUR 22-009-7,1.7. Man-Machine 
fntegration Branch, NASA-Anies Hescareh Center. 



aP-j 

,3Py 

P&(X)j 

iix 

a.\ 


Note that Ihe matrix Ai is the negative of the Hessian matrix (4) 
with the last set of terms omitted. An advantage of the method 
of scoring is that the matrix Aj can be calculated from the 
gradient itself, and the second partial derivatives of Py are not 
needed. When inverted, as it must be for the iterations, the 
inversi? becomes the Cranicr-Rao lower bound covariance matrix 
of estimation errors in the parameter ,x. See Curry (Note I ) for a 
more detailed discussion of these point.s. 
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MUNOML: The Executive Search Program 

As can be seen from Equations (3). (5). and ( 6 ), the only 
quantities needed to carrs- out the search function are the m x n 
quantities Py, and the m X nX fi quantities aPjj/axij. We have 
written MUN'OML to catty out the search operation using only 
the numerical values of the probabilities and the partial 
derivatives. These numerical values are provided by a user 
supplied subroutine which is called once in each iteration cycle. 
This provides an e.vtremely flexible tool for estimating 
parameters of grouped data. i.e.. for multinomial distributions, 
Thus, it will accommodate any behavioral model in which the 
stimuli and responses are drawn from discrete element sets as 
long as dP|j/3xk is continuous. The theoretical form of the 
underlying mstriburiom; within the model are contained in the 
subroutine. The range of distributions and models is almost 
limitless; One may treat continuous sensory models (signal 
detection theory); Thurstonian scaling; discrete sensory models 
(rectUinear ROC curves); piedictions of judgments and choice by 
the method of successive categories (more Thurstonian scaling). 
The distributions may belong to the same location/scale family 
(Gaussian, logistic, arc sine, etc.) or they may involve imbedded 
parameters (e.a., the beta distribution). 

In addition to the ma.ximization of the likelihood function 
with reject to the parameters, we have included a 
gcodness-of-fit test for the final param*iter estimates, given by 


x’ 



(ra - 
riPa 


, where rj = 




j=i 


df = m(n - 1) ■- E. 


(7) 


The degrees of freedom are obtained by finding the number of 
independent measurements and subtracting the number of 
parameters being estimated. We also have included another 
computation in which ail cells for which the c.vpected number of 
responses is less than four are pooled with adjacent cells and the 
degrees of freedom are reduced accordingly. 

MODELS IN SIGNAL DETECTION AND RECOGNITION 

Psychotogical Continuum (Thtustoitian) Models 

Rating seales-The method of successive categories: In the 
method of successive categories 'Bock ^ Jones. 1968), or the 
rating method in ps>chophysics (Green & S vets, 1966), the 
observed measures are Pij, the probability of eliciting response 
Rj, given stimulus Sj. These sensory-continuum models assume 
that each presentation of a stimulus maps to a point on an 
inietnal continuum, and these same probabilities, as predicted by 
the models, take the form 


Py = P(Rj/Sj) = P(Cj ^ U < Cj.f.x/Si) 

<7i 01 Ol 

“ F(Ztj+x) - F(Zij) 

i= 1 , 2 , . . „ m j = 1 , 2 , , , ,, n. ( 8 ) 

F(*) is the theoretical distribution function; and the unit 
deviates are 

= -ai - (9) 


Here is the jth criterion, level: Oj and .^<1 are. the mean and 
standam deviation of, the discriminal dispersion of the ith 




stimulus: and bi= lla^ and ai = ui/aj are an alternative set of 
parameters inuoduced for convenience. 

To remove the ambiguity of the origin and scale factor of the 
distributions, we assign 

b, = 1 . which fixes scale 

- ... - ( 10 ) 
a, = 0 . which ti.xes ongtn. 

Equations (10) are particularly useful when using rating scale 
methods in psychophysics, since one can arbitrarily assign the 
noise process to stimulus S, . and all distributions arc scaled with 
noise standard deviation so that bj becomes the ratio of noise 
standard deviation to signal .standard dei-jation. in the 
equal-variance case, aj is the normaiiacd separation between 
‘’naisc" and “signal” distributions, often called dj. 

.Muilqile S.NR pairs with equal variance: One of the new 
models that we have investigated is for decision behavior at 
multiple signal-to-noisc ratios (SNR) in which the e.xperitttental 
conditions were such that an equal-variance Gaussian model is 
apropos. Rather than work with criterion cut-off points, a it is 
mote meaningful to deal with the tog likelihood ratio which, for 
equal-variance Gaussian distributions, can be shown to be 

A(u) = di(bi g - >/itaj + ai+i)] i = 1 .3. , . .. rn - 1 ■ (1!) 

In the above e.xpression. we have used the convention that tfie 
stimuli are ordered by equal SN'R pairs, i.e., S, and S. ha'e the 
same SNR, S, and S 4 are at another SNR, etc. \Vhen this 
e.xpression for the log likelihood ratio is substituted into 
Equation (9) for the dcsiates with u = eg. we obtain 

Zi+l.J = j/d'i + d'i/2 i = 1.3.. ...m - 1 

fI2) 

.j - ^i.j/dj — d|/2 I = 2.3, , . .. n . 


If one assumes a constant likeliliood decision rule, then \g = 
for all i. 

Discrete Sensory Models 

The discrete sensory models lead to the prediction of 
rectilinear ROC curves (e.g., see Kraniz, 1969 & Luce, 
l963).These may be e.xpressed by 

N N 

Ph = Po + 'i^PiSCc - i) PfA = qo + .s AqjS(c - i) 

1=1 


(13) 


with the constraint 


N N 


where N is the number of limbs in the RtXi; Ph and PpA the 
theoretical values for the probability of a hit and false alarm ; p,, , 
Api, qo, Aqi describe the form of the ROC curve; e is the 
criterion level: and S(- ) is a saturation function. 


S(X)’ 


0 X ■« 0 

{x 0 < X < 1 , 

1 1 < .X 


04) 


Equations (13) and (14) arc suitable for use with the iiiethqd of 
scoring as an iterative technique since the gradient i> reasonably 
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v\i;ll bi;liuvcd. allhungh only picccwisc com inuouii with respect 
ui c. 

SUMMARY 

Uc liuvc found MUNOMl lo be cNtromeiy useful in 
da>*iutluy oporutiun since it provides the maximum amount i<i' 
llexib'tiity with a m'mtinum amount of duplicated elTort. This is 
especially true when minor changes to a model's structure are 
made because there veill be a great deal of unchanged code 
within the use) supplied subiouiinc, 

MUNOML has excetleni conveigence qualities and rarely goes 
beyond 10 ileraiionst when it does, it is usually because of 
programming errors in tl)c ealeulatiun of the piobabUiiies or 
their derivatives. 

VU are collecting a library of subroutines for use with 
.MUNOML. There are eutremly subroutines for calculating 
parameters for a successive categories model und for models 
based on the tlieury of signal detectability: unequal variance 
KU(' parameters from binary responses: equal-variance 
signal-to-floise ratio pairs; constant likelihood ratio decision 
rules; and Neyman-Pearson decisitm rules (usin;.’ objective and 
subjective probabilities). Cumpleie listings and/or curd decks 
t about SOD cards in all) are available on request. Operational 
details are described in Curry (Note 1). 

REFERENCE NOTE 

1. Curry, R. E. MUTCOML: A multinomial maximum 

likellbood proeram for behavioral research. Man-Vehicle 


L^oratoiy, Massachusetts Institute of TecbnoloBy, Caiubridee, 
Massachusetts, 1974. 
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COMPUTER TECHNOLOGY 

A random search algorithm 
for laboratory computers 


RENWICK E. CURRY 

Massachusetts Institute of Technology, Cambridge, Massachusetts 02139 

The small laboratory computer is Meal for experimental control and data acquisition. Postexperimental 
data processing is many times performed on large computers because of the availability of sophisticated 
programs, but costs and data compatibility are negative factors. Parameter optimization, which subsumes 
curve fitting, model fitting, parameter estimation, least squares, etc., can be accomplished on the small 
computer and offers ease of programming, data compatability, and low cost, as attractive features. A 
previously proposed random search algorithm {“random creep") was found to be very slow in convergence. 

We present a new method (the "random leap" algorithm) which starts in a global search mode and 
automatically adjusts step size to speed convergence. A FORTRAN executive program for the random leap 
algorithm is presented which calls a user supplied function subroutine. An example of a function subroutine 
is given which calculates maximum likeUhood estimates of receiver operating characteristic parameters 
from binary response data. Other applications in parameter estimation, generalized least squares, and 
matrix inversion are discussed. 

For many investigators involved in behavioral research, several of these are available in coded form (e.g., IBM 
the small laboratory computer is viewed primarily as an Scientific Subroutine Package), These algorithms 
experiment monitor, controller, and data acquisition typically require gradient calculations, some require 
device. The small computer is ideal for this purpose matrix inversions, and all would quickly overwhelm the 
because of its low initial cost and because it is dedicated capabilities of a small computer. As with most 
to the laboratory in which it resides. algorithms, they only attempt fo find one local 

The decision is not as clear when considering whether extremum, 
or not to perform postexperimental data processing on Other algorithms utilize a direct search process and 
the small computer. In favor of the small computer are ere more readily implemented on a small computer, 
data compatibility (no tape-to-tape or tape-to-card They are conceptually tt s most simple: the objective 
conversions), low operating cost, and availability, function or cost function f(') is evaluated at a point in 
Furthermore, there are many programs designed for parameter space a trial point 3£ = x + Ax is selected, 
these applications (e.g., the DECUS library). and if f(X) is an improvement over f(^), it is replaced by 

The advantages of processing the experimental data 5t. Thus, 'the best solution is always retained and 
on a large computer are the increases in flexibility, relinquished only when a better one is found, 
scope, and sophistication of processing techniques Algorithms of the direct search type are differentiated 
attendant with the increase in memory and by the manner in which the trial vMue (or Ax) is chosen, 
computational speed. There are many well-documented Examples of direct search algorithms are Chandler’s 
programs in existence which take advantage of these (Note 1) STEPIT and Hooke and Jeeves’ procedure 
attributes (e,g., the UCLA BMD series). (1961); these and other methods are compared in 

A strong case can be made for the small computer in Dorfman, Beavers, & Saslow (1973). The performance of 
one area of “sophisticated” data processing, and it has these algorithms, as measured by the conventional 
received only a modest amount of attention. Parameter yardsticks of computing time or the number of times the 
optimization, which subsumes curve fitting, model cost function is evaluated, is usually (but not always) 
fitting, parameter estimation, least squares, polynomial inferior to the more sophisticated approaches. For the 
root finding, etc., is a task which can be performed on small computer, however, time is not money, and when 
the large computer as well as the small. Many algorithms the further advantages of less programming effort (e.g., 
have been proposed (e.g., Fletcher & Powell, 1963) and no gradients to be computed) and data compatibility are 

considered, the direct search algorithms on a small 
This research was sponsored b> NASA Grant NGR 22-(X)9.733, computer become very attractive. 

Man-Machinc Integration Branch, NASA-AMES Research Center. One final feature, and some would say that it Is the 
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most important, is tiiat these techniques can more easily 
search for global, not just local extrema, since the logic 
for global search is readily incorporated in the direct 
search lo^c. 

In this paper, we present a direct search program for 
function minimization wliich is intended especially for 
use on small computers. Because of this goal, the 
primary consideration is small program size rather than 
speed of computation, number of function evaluations, 
or computation costs. The program’s genesis is the 
“random creep” algorithm originally proposed by 
Faureau and Pranks (1958) and discussed most recently 
by Bekey and Ung (1974), After describing the 
algorithm, we propose a new algorithm (the “random 
leap”) which offers substantially better convergence 
speed. A FORTRAN program to cariy out the random 
leap procedure is given, and an example of maximum 
likelihood estimation of receiver operating characteristic 
parameters is presented. Applications to generalized least 
squares and matrix inversions are also discussed. 

RAISTDOM SEARCH ALGORITHMS 
FOR GLOBAL EXTREMA 

The Random Creep Method 

The random creep algorithm (Bekey & Ung, 1974; 
Faureau & Franks, 1958) derives its name from the 
manner in which the trial values of the parameter vector 
increment Ax are derived: (l)At each stage of the 
iteration, each of the n elements of the currently 
optimal parameter vector it is perturbed by a zero-mean 
Gaussian random number of specified standard 
deviation: 

Ax| = e| 

j = 1 j . . . ,n 

ej =N(0,of) 

where Ax| is the perturbation to the jth element of x on. 
the ith trial, and the ej are Independent zero-mean 
Gaussian variables. (2) If the trial vector X - ^ + Ax 
results in an improvement, $ is replaced by X. If not, a 
new Ax is chosen according to the above equation. (3) If 
the iteration fails to improve after a specified number of 
consecutive trials, all standard deviations are increased 
by the same factor 

ffj=OiOj j=l,. ,.,n 

where ct> 1. This action is based on the assumption that 
a local miniirium has been reached. (4) The process 
terminates after either (a) the total number of iterations 
reaches a predetermined value, or (b) the standard 
deviations have been increased a specified number of 
times. 

The major advantage of this approach in small 
computer applications is the relatively small amount of 
storage space required for the search algorithm. There is 




Figure 1. Gcometiy of possible improvement with randpm 
steps (Ax) of flxed size (R) in two^parameter space, (a) Far from 
the extremum: the gradient vectors ^ are parallel, and an 
improvement is possible only if 0 (the ai^e between Ax and g) is 
-jr/2 < 0 < w/2, which occurs with the probability .5. (b) Close 
to a spherical extremum (E) at distance Pq : contours of constant 
criterion f(x) are concentric circles centered at E. Improvement 
is possible only if -0(, < 0 < 0^ < w/2 with probability < .5. 
When the step size is greater than twice the distance to the 
extremum (R > 2p,), no impiovement is posable. 

no guarantee that this algorithm will converge to the 
global extremum (only an exhaustive search will do 
that), but experience has shown this to be an effective 
method of finding local and global extrema (Bekey & 
Ung, 1974). 

The Random Leap Method 

Improving convetgence rate. The main disadvantage of 
the random creep method, as it is described In the 
previous section, is the relatively slow rate of 
convergence because only the minimum step size is used. 
We experimented with modifications to rectify this 
problem; for example, if Ax had been an improvement 
on the previous trial, it Was then used in succeeding trials 
until no further improvement was obttoiied. These 
modifications gave a minimal amount of improvement. 

It was about this time that the analytical work of 
Rastrigan (1963) came to our attention; He analyzed a 
random search procedure where the direction of the 
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random search vector is uniformly distributed over the 
hypcrspliere, but the step size is fixed. The areas of 
interest are shown in Figure I fora two-parameter case; 
Ax is the random increment in the parameter vector and 
has a fixed magnitude or step size of R; 0 is the angle 
between the gradient vector g and the random increment 
Ax. In Figure la, the currently optima! value of x is far 
from the extremum resulting in (nearly) parallel 
gradients througlio’ut the trial region. The probability of 
an improvement is under these conditions. 

Figure lb shows the situation when the search 
procedure nears an extremum E which is assumed to be 
of a spherical nature-the contours of constant cost 
function are circles in this two-parameter case, but are 
hyperspheres in general. In this diagram, the step size is 
again denoted by R; 0 is the angle of Ax from the 
gradient, 0o is the maximum angle for which an 
improvement can be made; and pQ is the distance of the 
current point from the extremum. It can be seen that 
the probability of improvement is less than H under this 
condition, since there is a smaller region in which the 
random search vector will yield an improvement. In the 
extreme case, when the step size is twice the distance to 
the extremum, then all points of the function which 
would yield an improvement lie within the circle of the 
search vector and no improvement is possible, These 
factors are shown quantiatively in Figure 2, which 



8TEPS1IE <R)/DISUItC8 FIW EITBEMM (p^> 

Figure 2. FtobabiUty of no improvement vs. step size 
(measured in distance to extremum) for spherical extremum in 
n-dimendonal space. 
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F^re 3, Approximate magnitude R of the random search 
increment ^x. Upper abscissa scale for R normalized by the 
standard deviation (Gaussm distribution); tower abscissa scale 
for R normalized by the half-limit of a uniform distribution 
(DLM) = \/3a). Circles show exact values for uniform 

distribution (n ~ 2). 

demonstrates the probability of no improvement for 
fixed step size as a function of the step size, measured in 
distance from a spherical extremum, for various 
dimensions of the parameter vector. These curves were 
generated by numerically integrating the expressions 
given in Rastrigan (1963). 

Figure 2 contains the information that is the basis for 
the random leap method. In essence, the step size is 
adjusted so that one always obtains a modest probability 
for improvement. When the step size is large, it is 
difficult to get an improvement, so the step size is 
halved. As the search procedure then comes closer to the 
extremum, it again becomes difficult to obtain an 
improvement and the step size is again halved. 

When dealing, as we are, witli a completely random 
increment in each coordinate direction, the step size is 
not fixed. In Figure 3, we show the cumulative 
distribution curves for the distribution which is the 
magnitude of the random increment Ax when the 
increments are independent Gaussian variables with the 
same standard deviation. The algorithm we suggest uses 
uniformly distributed variables, so we have also shown, 
for later use, the abscissa in the scale of step size 
normalized with respect to the limits of the uniform 
distribution. The difference between the uniform and 
Gaussian distributions (in terms of the magnitude of the 
vector) becomes small for a large number of parameters; 
the circles in Figure 3 indicate the exact solutions for 
uniformly distributed variables for n = 2, and even for 
this low number of parameters, the deviations are of no 
practical significance. 

Program outline. The overview of the operation begins 
with a call to the user-supplied function routine to 
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CrrjcuMAAKOQH leap pAftAMEreq SEAACH ftaLT|h<£ 

04KENS10M X(5?.TMPXC5;,PLH0(>),TWLH«5> 
fai:sBuj£T J At:«MCPy 23- BNDMfi=J«“<P-l)" I FOft P EIT MACHINE 
lA s 6l 
RNPHQ a 2g(i7 

IhND = 1237 

taMMHHiHjTULUATlON CALL TO CET PAflAMEfCRS ANO PlHST VALUE 
1 LfR$T s I 

TCST=CSTFNCTHPX,NX,LFRST,H>CJ TR,MXFGS*MXEXP,MXFLS,J1HCTA,LP(IT^PLM03 
LFR$T ^ 2 
NPAIL = 0 
I TEA = 0 

CSB53CSTAAT LOCAL SEARCH HERE 
SO HOPE « t 

C“=w»«START global search here fWOOE SET TO 2 SELOW) 

33 NEXP 3 0 
NCTR a 0 

CMKBKUsET TEMPORARY OJSTfi. LIMITS AT tHJTJAL VALUE 
DO Ld I si, NX 
‘♦fl TMPLMCI) = OLMOCl) 

GO TO CSQ/6Q),H00S 

C“M««»^SAVE SUCCESSFUL TRY AND PRINT ]F LPRNTs? 

50 CST = Tcsr 
OD 60 i=l>NX 
6fl XCl) = THPXCJJ 

CO TO C75,703jLPRT • 

70 UftlTE |TER,KODE,NFAJL.NE<P,WCTR,CST,Cx<13,| = 1,NX) * 

1Q40 FORMAT C 3J 6, E 15 . 3, EFIO .b/ C lOX, Id > 

75 NFAtL = 0 

Cusssapjup Xejl TRIAL X AFTER ChECXINR FOR MAX ITERATIONS 
SO ITER £ ITER t 

If O'iCa - MXiTRi 90,^0, 3Jfl 
90 DO 99 lsi,NX 
IRNO = lA^IRND 
IF <IRND) 99,230,99 

99 THPXCt) = XCI) FLOAT(lRND>/RNDMjJ‘TMPLHCl> 

TCST3C5TFNCTMPX,MX,LFRST,NxI T ft,MXfCS,MXEXPjNXFLS,HXCrK,LPflT,OLMD) 
IF CTCST - CST> 100,510,110 

CSSSiisrsAVC IMPROVED SOLUTION OUT FIRST SWITCH TO LOCAL IF |N GLOBAL 
loo GO TO C30,30J,MDDE 

C»“«»»CTHEAYII$E COUNT SUCCESSIVE FAILURES, 60 TO STEPSfZE LOGK 
110 NFAIL £ NFAIL f 1 
CD TO O20, 1703, MODE 
C««=:«W5T£PS[ZE LOGIC FDR LOCAL SEARCH 
Ctt««aWTE5T FOR MAX FAILURES, LOCAL SEARCH 
120 IF {NFAtL . MXFLS) 90,9D,iSD 
130 NCTR s NCre 4 1 

IF CnCTR - MXCTR3 ILQ, Itifl, 160 
Cwnr.aSHALVE DISTRIBUTION LIMITS AND TRT NEW X 
ibO NFAIL £ 0 

00 ISO |sl,NX 
130 THPLrtCl) = TMPLMCn-.S 
GO to 9D 

CHWHSIBeND local search if too MANY CONTRACriCNS CMALV]NCS> 

160 WRITE a, 1100) 

HOC FORMAT C^gNO LOCAL SEARCh'7 

WRITE Cj,lgbO) tTEfl,M0De,NFAIL,hexP,t<CTa,CST,CXCl3, l = l,NX> 

Mooe - 3 
CD TO 35 

C«a»Jf«5TEPSIZE LOGIC FOR GLOBAL SEARCH 
Ca«:«:iJiTEST FOR MAX FAILURES- CLCflAL SEARCH 
170 IF CNFAIL « HXFGS) 80,90, iSO 
tao NEXP a NEXP •» 1 

IF CNEXP - MXEXP3 190,190,3)0 
Caaa;jS£j,pAHO DISTfitBuriON LJHITS 
190 NFAIL £ 0 

DO 200 }Kl,NX 

200 ThPLM(I) s 1, J-’TMPLmCO 
GO TO 80 

CuaaHUgyiT SEARCH LOGIC HERE AND START AGAIN 
730 WRITE a,lDB0) 

1060 FORMAT {'STOP') 

WRITE {l,10ba) iter, MODE, NFAil,NEXP,NCTR, CST, C x{0,j£l, NX) 

GO to I 
END 

Figure 4 


supply the constants for the program operation and the 
initial guess at- the parameter vector. From there, the 
operation proceeds as follows: (a) Uniformly distributed 
independent increments of each coordinate are chosen 
to calculate the new trial value of x.^ Initially, the 
distribution limits of these random increments are 
chosen to be of moderate size relative to the entire 
search region, Le., if the search region is (—1,1) in each 
coordinate, then the limits of tlie uniform increments 
are on the order of (—1,1). Tliis allows a preliminary 
global search in the beginning phases, (b) When the trial 
values of x have failed to yield an improvement a 
specified number of times, it is assumed that the step 
size is too large relative to the distance to the extremum, 
and the distribution limits are halved. The search 
continues with these limits until the procedure fails to 
yield an improvement in the (same) consecutive number 
of trials, (c) After the distribution limits have been 


lialved a specified number of times, it is assumed that 
tfie procedure has converged to a local minimum. The 
results are printed out, and tiie distribution limits are 
reset to their original, moderately sized, values. The 
global search is then initiated in a manner similar to the 
random creep method: if no improvement is reached 
after a specified number of times, the distribution limits 
are increased, and the search continued. Note, liowever, 
that the glob.d search is initiated witli tlie moderately 
sized distribution limits and not the minimum size, as is 
done in the random creep melliod. (d) The procedure is 
terminated whenever the total number of iterations 
exceeds a specified value or when the global search lias 
not yielded an improvement after a specified number of 
expansions of the distribution limits. 

The advantage of the random leap algorithm lies in its 
ability to perform a preliminary global search and 
gradually reduce the step size as tiie extremum is neared. 
There are cases, obviously, when the random creep 
method would be better, such as wlien tiie initial value 
of x is very close to die extremum. 

Program description. A FORTRAN listing of the 
random leap algorithm is given in Figure 4, wiiich 
executes the general procedure described above. To do 
this, it requires a user-supplied function subroutine 
CSTFN wliich calculates the value of the cost function; 
we have made provision for the user to enter the 
constants and parameters for the search routine for 
maximum flexibility. During the first call to the 
function routine CSTFN, the following transfer of 
variables obtains: 

Variables Passed to CSTFN 

LFRST has been set to the value of 1 to indicate Uiat 
tills is the first (initialization) call 

Variables Returned from CSTFN 
TMPX temporary value of X, in this case the initial 
value 

NX the number of dimensions in die parameter 
vector X 

MXITR maximum number of iterations 
MXFGS maximum number of consecutive failures in 
the global search mode 

MXEXP maximum number of expansions (increases of 
the distribution limits) in the global search 
mode 

MXFLS maximum number of consecutive failures in 
the local search mode 

MXCTR maximum number of contractions (halvings) of 
the distribution limits in the local search mode 
LPRT returned as 2 to print ''a) each time an 
improvement is made, (b) at the end of the 
local search, and (c) at the termination of the 
search; returned as I to print only at the end 
of the local search and at die end of the entire 
search 
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TCST ilie numerical value of the cost function 
DLM0 the array containing the initial limits of the 
uniformly distributed random search 
increments. For example, AX(l) is uniformly 
distributed between (-DLM0(1), DLM0(I)), 
For subsequent calls to the function routine during 
tlie normal course of operation, the following transfer of 
variables applies; 

Variables Passed to CSTFN 
TMPX the temporary or trial value of X 
LFRST has been set to 2 to indicate that initialization 
is not required 


Search, sfiould be in the range of 9 to 13, perhaps 
greater to account for contingencies. 

The parameters for the global process (MXFGS and 
MXEXP) can really only be determined by 
experimentation in each particular situation. We found, 
in the cost function described below (which has two 
maxima) that 300 trials with a five-dimensional vector 
were more than adequate to find the global maximum 
when situated at the local maximum. The number of 
such trials must be increased substantially when going to 
a larger paiameter vector to obtain similar densities of 
trial values within the parameter space. 

APPLICATIONS 


Variables Returned from CSTFN 

TCST the numerical value of the cost function 

CHOOSING THE SEARCH PARAMETERS 

In this section, we give some guidebnes on how 
one determines the parameters of the search. Assume 
that the parameters have been scaled such tliat the 
solution is likely to lie within the unit hypercube (-1,1) 
in each coordinate. Let us pick a convergence criterion 
that (say) eadi value sliould be determbed withm at 
least .01 of the true value which mmimizes the cost 
function. To find the minimum step size required, we 
arbitrarily fix the probability of no improvement which 
we would like to detect at (say) .7. Looking at Figure 2, 
we can then determine the step size relative to the 
distance from the extremum for various dbiensions of 
the parameter vector. If we have four parameters, then a 
probability of .7 occurs with a step size that is roughly 
.6 times the distance from the extremum. Smee we want 
the distance from the extremum to be no more tlian .01 , 
the magnitude of the minimum step size must be .006. 
To find the mbimum distribution limits corresponding 
to the mbimum step size, we refer to Figure 3, where 
we see that for n = 4, the step size will be less than 1.8 
times the distribution limit more than 95% of the time. 
Thus, the mbimum limits of the uniform distribution 
should be ±,006/1.8 = ±.0033. The bitial value of the 
distribution limits should be on the order of ±1, and 
because tlie limits are halved at each stage, we seek an 
integer M such that 1/2'^ .0033. Note that 

2®(.0033) = .768, so we set DLM0 = .8 and 
MXCTR = 8. 

The adaptive nature of the random leap 

program results from testing the hypothesis that the 
probability of failure is .7 or less. To reject this 
hypothesis at the .05 and .01 levels of significance, we 
would need 9 amd 13 successive failures, respectively, 
(six and nbe successive failures are requited for the same 
level of confidence when detecting probability of failure 
,6 or less, showbg the desirability of workbg with 
probabilities of no improvement closer to .5 than 1 .0.) 
Thus, MXFLS, Maximum Consecutive Failures, Lpcid 


Maximum LiltelQiood Estimation 

In tliis section, we describe an application of the 
random leap algorithm and the results of a Monte Carlo 
simulation. The first example uses a cost function of tlie 
form 

f(x) = .59exp(-Qi (x)) ± ,41exp(-Qa(x)) (1) 

Q, = S(n) |^(.5 - /(1 .26)‘-i (2a) 

Qa = S(n) |^(.5 + xO® /(1 .26)t-i (2b) 

where S(n) is a scale factor dependbg only on the 
number of dimensions b the parameter vector and keeps 
tlie argument of the exponent function williin 
reasonable bounds. This function was chosen for a 
detailed exambation because of its two nearly equal 
modes. It corresponds to a maximum likelihood 
estimation problem ta which the observations (.5, 
.5, , , .) come from either a distribution with mean x 
(prior probability .59) or — x (prior probability .41) and 
unequ^ variances. 

Monte Carlo trials. One hundred Monte Carlo trials of 
the random creep and random leap operations were run 
usmg the cost function m Equation 1 on a PDP-12 
computer with software multiply and divide routines. 
Approximate time for each iteration was about l.S sec. 
The bitial conditions for these trials were uniformly 
distributed in the unit hypercube, and parameter vectors 
of dimensions 2 and 5 were used. Convergence was 
considered complete when the search routine first came 
withm a sphere corresponding to an rms deviation in 
each component of .05. 

With the two-dimensional parameter vector, the mean 
number of iterations-to-convergence for the random leap 
algorithm was 35, the median 36, and the maximum 
number of iterations required was 67. In all these cases, 
the fact that the search was started at the global level 
rather than the local level resulted in a convergence to 
the global maximum first. On the other hand, the 
random creep algorithm converged to the local 
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maximum first in many trials. It also failed to find the 
global maximum within the allowable maximum of 500 
iterations on many of its trials, whereas the random leap 
algorithm never failed to converge to the global 
maximum. 

In testing the random leap algorithm with the 
five-dimensional parameter vector, the mean number of 
iterations-to-convergence in 100 Monte Carlo trials was 
194, with a maximum number of iterations of 365 and a 
minimum of 58, The distribution of the number of 
iterations had two modes; one, at 100 iterations, 
corresponded to those searches that went directly to the 
global maximum; the other mode, at 240 iterations, 
represented those searches that went to the local 
maximum first and then to the global maximum. Sixty 
of the 100 Monte Carlo trials converged to the global 
maximum first; this is a statistically significant 
difference from the expected value of 50 tliat one would 
expect with tlte random creep algorithi.i, and is due to 
the fact that the random leap algorithm starts off in the 
global search mode. The random creep algorithm was 
not run in this more difficult task because of its 
relatively poor performance on the easier two-parameter 
case. 

Receiver Operating Characteristics 

In the theory of signal detection (TSD) approach to 
psychophysics, the subject’s response is divided into 
sensory and response bias components. One may 
administer a yes-no response procedure in several 
sessions, during each one of which a subject adopts a 
different strategy or criterion level (|3), whereas the 
sensitivity of the signal {d') remains constant throughout 
(Green & Swets, 1966), Assuming that we have binary 
response data (“YES“ or “NO”) for M criterion levels, the 
TSD model uses the following expression for the 
probability of a “YES” response given noise (n) and 
signals) respectively 

P(“YES”lnoise) = 1 - ^(0i) = ^(-ft) (3) 

i = l,2, ...,M 

P(“YES”lsignal)= 1 -^[(0i -d')/Osl =‘b[(d'-/3,)/f7j 


where Pi is the criterion level adopted in the ith session, 
d' is the sensitivity measured in noise standard deviation 
units, Oj is the standard deviation of the signal 
distribution, measured in noise units and 4* is tlie 
Gaussian distribution function. 

The maximum likelihood estimates for the parameters 
01 , d', and Og are obtained by maximizing the likelihood 
function which is proportional to 


m 

C(x) « f(x) = r,„„In|cP(iSjj + r,YnInj4’f-0i)l 


+ r,„glnjtl»(/3i - d')/(7g[ + rjYsln{4>(d' -/3j) /Ogj (5) 


where rij^ is the number of j responses at criterion i to 
stimulus k, i = 1,2, ... Al;j = YES,NO;k = s (Signal), n 
(noise). The parameter vector for this case is 



( 6 ) 


Note that we can easily obtain the parameter estimates 
under the constraints of equal-variance distributions by 
setting Og in the above expression to unity. 

A function subroutine written to generate the 
parameter estimates is shown in Figure 5. It prompts the 
user for the number of different criterion levels, and a 
parameter indicating whether it is to be an 
equal-variance case or whether B = I/Og is a parameter to 
be estimated. ^ 


PUNCTIOH CSTPNCX^NX^LPAST^HXt TR^HXF MXEXP, HXPLS, MXC TA, LPRT^ 0 LHQ > 
^:mu»upij^CT10h &UBAOUT1NE POA MAXIMUM LlKEL^ .OOD ESt|HATlOM QF RQC PARAH 
CMitiritifKJitTTEN FOR TELETYPE INTERACTION 
DIMENSION XCl),AC3,fi),OLMO< O 
CO TO Cl,30),LFA5T 

CB«nnit7HIS PORTION FOR TmE INITIALIZATION CALL TO INPUT OATA 
1 NAITE CUIOOO) 

lOaO FORMAT O RE'w 'dX^HF;TR,MXFGS,MXEXP,HXFLS,MXCTR,LPAT*5 
READ NX„nx| TR^MXPOS^M xEXP^MXPLS^MXCTR^LPRT 

1010 FORMAT 

WRITE loos) 

lOOS FORMAT C NUMBER OF CRITERIA/! FOR EQUAL VAA^T OTmERnICE*) 

READ Cl^tOlO) HCAIT.LO 
IF 1,5.1 

5 WRITE a, lOTO) 

lOTO FORMAT READ IN0,«YES RE5PDNEES;FIR5T NOISE, THEN S *■ N') 

DO 10 tCAIT ? 1,MCRIT 
OD U tsic » 1,3 
I « 1519 f Zi^ClCRIT - O 
10 READ Cl, 1030) AC 1, 1 ) , RC3, 1 ) 

1030 FORMAT (F)0.2) 

WRITE Cl, 1040) 

Ui*iQ FORMAT V READ XOCl )/ DLHOC 1 3 ’ > 

IMX & HCrIT ^ LQ 
OD 30 I 3 1,IMX 
30 READ C1,IQ5Q> XC{),DUH0CI) 

CWW^-WWnOAMAL calculations kerb 
C» wnwU5TRIP D* FROM X AND LIMIT IF NECESSARY 
50 IX 3 MCRIT 4> I 

IF Ci^Clx) > .01) 34,56,58 
5A XCIX) 3 .01 
36 0 B XCJX) 

C»hk«»5bT STD. DRV. RATIO 
CO TO C*»0,50),LB 
40 Q B 1. 

CD TO 6l 

50 rx 3 HCRIT * 2 
C«h»«hi.|MIT ALLOWA&LB Q ,CE. 0,01 
IF CXCIX) - .01) 54,58,58 
34 /CIX) = .01 
SO e 3 xcix) 

C»fiRAHCALCULATE LOO LIKELIHOOD FUNCTION <LE55 CONSTANTS) 

60 XILP s 0. 

DO lOD tCRlT = 1,MCR|T 
OD tac t$io B i,a 
CO TO C7O,80),rslC 
70 2 5 XCICRIT) 

CO TO 80 

80 2 B e«CxCtCRlT) - D) 

CXUNWFTKEOfleTtCAL PRD6C*NO'> PROM CAUS5IAN CUMULATIVE DIST.TIBUTION 
90 P CCOFCZ) 

IP = ESIC * 2WCICRIT 1) 

100 XLLP 3 XLLP 4- RCl, IP)«AL0CCP) RC2, |P)»ALOCO . • P) 

C«»*«wminIH1ZE the negative LOG LIKELIHOOD FUNCTION 
CSTPH 3 • XLLP 

return 

END 

FUNCTION GCDFC7) 

AZ s APSCZ) 

T = t./Cl. 4- ,aSl6419WA2) 

D : *59S9425»EXPC - Z»Z».5) 

P B 1. - 0«T«CCCC1.33fl374«T - 1,0J1256)«T ♦ 1,781478)«T 
Q - .356SiS8)»T + ,3193815) 

IF CZ) 1,2,2 

1 P 2 1, « P 

2 GCOP = 9 
RETURN 
END 


Figure S 
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For a test case, we assumed a value of x equal to 
(Pt = .5, ^2 = .75, 03 = 1.25, d' = L00,(7s=i-25). We 
used as responses the values proportional to the 
theoretical probabilities, i.e., we expect a “perfect” Fit. 

^ Starting from the initial condition (-.5, 0, .5, ,5, .5), the 
random leap algorithm converged to within .01 of the 
true value of each component of the parameter vector 
within 100 iterations. To determine how well these data 
from an unequal variance distribution are fit by an 
equal- variance model, we used the equal-variance option 
of the program and found the equal-variance fit to the 
unequal-variance data; the final estimates were 
(0, =.531, 02 =.756, 03 = 1.1, d' = .96), and the 
random leap algorithm, starting from the same initial 
conditions as before, again converged to within ,01 of 
the final values after 100 iterations. Each of these 
solutions took about 3 min on the PDP-12 computer 
described above. Since a printout was made at each 
improvement, though, the computations were limited at 
limes by the print operation. 

Parameter Estimates from Grouped Data 

The maximum likelihood procedure is one method of 
obtaining estimates of distribution parameters for 
grouped data, for example, obtaining the mean and 
standard deviation of a distribution from data such as is 
gathered in poststimulus histogram form. In these cases, 
the log likelihood function is proportional to 

E(x) oc LilnPi(x) (7) 


calculating the value which can be used in a 
goodness-of-fit test when the final estimates have been 
obtained. 

Generalized Least Squares 

One of the most common and most powerful forms of 
regression analysis is generalized least squares, which can 
be written in the form 

f(x) = (y - h(x)FW6^ - h(x)) (10) 

where y is a k-dimensiona! vector of observations (data), 
h is a vector function of the parameter x, and W is a 
positive semidefinite (symmetric) matrix of weighting 
coefficients. The objective is to choose the parameter 
vector X to minimize this cost function. In the linear 
form, Equation 10 becomes 

f(x) = (y - Hx)TW(y - Hx) (11) 

where H is a kxn matrix. This has a unique solution 
under suitable conditions on H and W (Rao, 1973), 
ususally satisfied in practice. One may use the random 
leap algorithm to find the solution to Equation 1 1 rather 
than go through the matrix inversions and manipulations 
required to solve for the vector x. We also note that if H 
is square and nonsingular, and if we solve Equation 11 
with y equal to zero except for a 1 in row i, then 
solution x is the ith column of the matrix H" L Thus, n 
separate parameter searches will completely invert the 
nxn matrix H. 


where r, is the number of responses in the itli group 
(i = i,2, . . . ,N), and Pi is the theoretical probability of 
the sample falling in the ith group. Under the 
assumption that the underlying distribution is Gaussian, 
Pi is given by 

Pi = ^{(ci+ 1 - M)/ff} - *|(ci - (i)lo] (8) 

Oi = cn+i= 

and the parameter vector is 

x=(g) (9) 

Simple changes can be made to the function routine 
shown in Figuie 5 to perform these calculations. The 
boundaries of the groups jcijmust, of course, be read in 
during the initialization phase, but it is a simple matter 
to program Equation 7 and solve for the mean, fx, and 
standard deviations, a. 

An alternative to using the maximum likelihood 
estimation criterion is the minimum criterion (Rao, 
1973). This has many of the advantages of the maximum 
likelihood procedure, and serendipitously, one is 


CONCLUSIONS 

The disadvantage of a small laboratory computer for 
postexperimental data processing is its small memory 
and inability to accommodate large sophisticated data 
processing programs, while the advantages of using the 
small computer are data compatibility, ease of program 
development, ability to run for long periods of time, low 
cost, and accessibility. Parameter optimization is one 
area which has received relatively little attention in small 
computer applications. Direct search methods are the 
easiest to program, and they require little core storage 
and no analytical gradient calculations. 

Random search algorithms are one member of this 
class, and Rastrigan (1963^ has shown that the average 
rate of convergence of the random creep algorithm may 
actually be superior to the gradient method in which the 
gradients are calculated numerically at each step. The 
slow rate of convergence experienced in the random 
creep algorithm is due to small step size and was 
overcome by the random leap algorithm proposed here. 
This algorithm operates by starting off in the global 
search mode and automatically reducing the step size as 
the search procedure approaches the extremum of the 
function. After reaching a minimum (or maximum), it 
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branches out and starts a global search by gradually 
expanding the size of the random increments from their 
original, moderately sized values. 

The random leap algorithm was compared to the 
random creep method on a bimodat likelihood function 
and showed superior convergence characteristics. We 
then presented a function subroutine to be used with the 
random leap program to calculate the maximum 
likelihood estimates of receiver operating characteristic 
parameters in YES-NO tasks; this was followed by a 
description of maximum likelihood estimation of 
distribution parameters from grouped data (histograms). 
Generalized least squares regression was also considered, 
and it was shewn how one could perform matrix 
inversion using successive applications of direct search 
methods. 
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repcuting the sequence. 
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